Summary:
Analysis log for clinal fundulus project

Analysis Outline

grViz("digraph flowchart {
      # node definitions with substituted label text
      node [fontname = Helvetica, shape = rectangle]        
      tab1 [label = '@@1']
      tab2 [label = '@@2']
      tab3 [label = '@@3']
      tab4 [label = '@@4']
      tab5 [label = '@@5']
      tab6 [label = '@@6']
      tab7 [label = '@@7']
      tab8 [label = '@@8']

      # edge definitions with the node IDs
      tab1 -> tab2 -> tab3 -> tab4 -> tab5 -> tab6 -> tab7
      tab6 -> tab8
      }

      [1]: 'library prep (Elshire GBS)'
      [2]: 'sequencing (illumina hiseq2000)'
      [3]: 'qc checks and datawrangling (fastqc and bash)'
      [4]: 'demultiplex and read filtering (stacks: process_radtags)'
      [5]: 'read mapping (bwa_mem)'
      [6]: 'genotype likelihoods (ANGSD)'
      [7]: 'demography and population structure'
      [8]: 'selection'
      
      ")

preserveb9b4482ed50a56cf

Data description and QC

sequence data

The clinal dataset is composed of 9 lanes of HiSeq2500 run in SE100. These reads include 1478 individuals from 5 separate experiments plus and additional set of lanes used for undersampled regions (CT, NC, SC and ME).

Lane names:

H7T6MADXX_1 H7T6MADXX_2 H9AREADXX_1 H9AREADXX_2 H9CWFADXX_1 H9WYGADXX_1 H9WYGADXX_2 HMMKLADXX_1 HMMKLADXX_2

Samples

sites<-read.table("./sites/sites.txt", header = T, sep = "\t")
kable(sites)
Population_Full_Name Abbreviation N Longitude Latitude
Brayton Point, MA BP 50 -71.18604 41.71250
Clinton, CT CT 24 -72.52448 41.27902
Sapelo Isalnd Ferry Dock, Brunswick, GA GA 41 -81.35452 31.44932
Horseneck Beach, MA HB 50 -71.02556 41.50449
Kings Creek, Severn, VA KC 25 -76.42462 37.29845
Mattapoisett, MA M 31 -70.81464 41.65875
Mount Desert Island Bio Lab MDIBL 24 -68.29441 44.42901
Wiscassett, ME ME 27 -69.66481 43.99695
Mantoloking, NJ MG 331 -74.07045 40.05046
New Bedford Harbor NBH, F, H, P, Syc 150 -70.90760 41.62486
Manteo, NC NC 24 -75.66737 35.90029
Oyster Creek, NJ OC 47 -74.18475 39.80866
Rutgers, Tuckerton, NJ RB 423 -74.32448 39.50878
Charleston, SC SC 26 -79.91624 32.75873
Stone Harbor SH 62 -74.76492 39.05666
Slocum River, MA SLC 30 -70.97109 41.52257
Succotash Marsh, Matunuck, RI SR 49 -71.52557 41.38235
Wilmington, NC (Carolina Beach State Park) WNC 30 -77.91932 34.04998
grandis grandis 34 -84.17740 30.07200

QC reports


#!/bin/bash

# set the number of nodes
#SBATCH --nodes=1

# set the number of cpus
#SBATCH --cpus-per-task=6

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-02:00:00

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of job
#SBATCH --job-name=fastqc

# set name of output file
#SBATCH --output=fastqc.out

# mail alert at start, end and abortion of execution
#SBATCH --mail-type=ALL

# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu

files="H7T6MADXX_1
H7T6MADXX_2
H9AREADXX_1
H9AREADXX_2
H9CWFADXX_1
H9WYGADXX_1
H9WYGADXX_2
HMMKLADXX_1
HMMKLADXX_2
"

for file in $files
do
/home/jgibbons/SOFTWARE/FastQC/fastqc -t 4 /home/ddayan/fundulus/seqdata/${file}_fastq.gz -o /home/ddayan/fundulus/seqdata/fastqc

done

fastqc looks good

Demultiplex and Read Filtering

We will use ANGSD to create population level allele frequency estimates rather than call indiviudal genotypes. The input files for ANGSD of indivudal level BAM files. Therefore the first step after QC check is to demultiplex, remove low quality reads and trim adapter and barcode sequences. Ths is accomplished with the process_radtags program from stacks v2.3.

generate barcodes

Barcode keys (files used to demultiplex sequence data) are a bit complex given that some fish are sequenced twice and stacks requires each lane be run separately through process radtags. The solution is to give duplicate fish a unique ID suffix, but with a shared prefix, then cleaned demultiplexed reads from each fish in each lane are combined with their additional reads using prefix matching.

make new fish_ids, based on only population and library prep id, then mark duplicates (duplicate fish_ids occur because some individual fish are sequenced twice i.e. in different lanes)

barcodes<-read.table("./snp_calling/barcode_keys/cline_barcode_key.txt", sep = "\t", header = T)
tmp<-paste(barcodes$Pop , "_", barcodes$LibraryPrepID, sep = "")
barcodes$fish_id<-make.names(tmp, unique = TRUE)
write.csv(barcodes, row.names = FALSE, file = "./snp_calling/barcode_keys/cline_barcode_key.csv", quote = FALSE)

used python script on hand for slightly different format: create a unique single variable for sequencing lane (i.e. colbind flowcell and lane with a _) then run python script below


""" The input file has four columns, this script takes writes columns 2 and 3 (barcode and individual) to a new file based on the value of column 4."""

import csv

with open('./snp_calling/barcode_keys/cline_barcode_key.csv') as fin:    
    csvin = csv.DictReader(fin)
    # Category -> open file lookup
    outputs = {}
    for row in csvin:
        cat = row['Flowcell_Lane']
        # Open a new file and write the header
        if cat not in outputs:
            fout = open('./snp_calling/barcode_keys/{}_key.csv'.format(cat), 'w')
            dw = csv.DictWriter(fout, fieldnames=csvin.fieldnames)
            dw.writeheader()
            outputs[cat] = fout, dw
        # Always write the row
        outputs[cat][1].writerow(row)
    # Close all the files
    for fout, _ in outputs.values():
        fout.close()

oops, only meant to keep barcode and sample id and need to write to tab delimited file


for i in ./snp_calling/barcode_keys/*csv
do
  cut -d "," -f 2,10 $i > ${i%.csv}.tmp
done


for i in ./snp_calling/barcode_keys/*tmp
do
    tr "," "\\t" < $i > ${i%.tmp}_barcodes.txt
done

clean up directory

rm ./snp_calling/barcode_keys/*.tmp
rm ./snp_calling/barcode_keys/*[12]_key.csv

remove first line from all files

sed -i -e "1d" ./meta_data/*_barcodes.txt

process_radtags

Here we clean and demulitplex the reads.

process radtags options:

-f single end
-c remove any read with an uncalled base
-q remove any read with low quality
-r rescue barcodes
–inline-null barcodes inline
-e aseI
–adapter-1 trim adapter sequence (used AGATCGGAAGAGCCGTTCAGCAGGAATGCCGAGACCGATCTCG (common adapter from elshire))

#!/bin/bash

# set the number of nodes
#SBATCH --nodes=1

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of output file
#SBATCH --output=library_%a_process_radtags.out

# mail alert at start, end and abortion of execution
#SBATCH --mail-type=ALL

# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu


source /opt/stacks-2.3/bin/source_me
/opt/stacks-2.3/bin/process_radtags -f ./seqdata/HMMKLADXX_2_fastq.gz -b ./meta_data/HMMKLADXX_2_key_barcodes.txt -o ./cleaned_tags -e aseI --inline-null -c -q -r --adapter-1 AGATCGGAAGAGCCGTTCAGCAGGAATGCCGAGACCGATCTCG  &> ./cleaned_tags/pr_library_HMMKLADXX_2.oe

Summary of process radtags:

A summary of number of reads for each lane is in the code chunk below

H7T6MADXX_1  
162739994 total sequences   
 13628228 reads contained adapter sequence (8.4%)  
   312933 barcode not found drops (0.2%)  
   828040 low quality read drops (0.5%)  
 14062227 RAD cutsite not found drops (8.6%)  
133908566 retained reads (82.3%)  
  
H7T6MADXX_2  
151505038 total sequences  
 13951706 reads contained adapter sequence (9.2%)  
   132166 barcode not found drops (0.1%)  
   584797 low quality read drops (0.4%)  
  8662501 RAD cutsite not found drops (5.7%)  
128173868 retained reads (84.6%)  
  
H9AREADXX_1_fastq.gz  
148867111 total sequences   
  4339944 reads contained adapter sequence (2.9%)  
   306024 barcode not found drops (0.2%)  
   305700 low quality read drops (0.2%)  
  5881380 RAD cutsite not found drops (4.0%)  
138034063 retained reads (92.7%)  
  
H9AREADXX_2_fastq.gz  
122597494 total sequences  
   720840 reads contained adapter sequence (0.6%)  
   185795 barcode not found drops (0.2%)  
   198200 low quality read drops (0.2%)  
  8394213 RAD cutsite not found drops (6.8%)  
113098446 retained reads (92.3%)  
  
H9CWFADXX_1  
110486132 total sequences  
 13263006 reads contained adapter sequence (12.0%)  
    98479 barcode not found drops (0.1%)  
   775153 low quality read drops (0.7%)  
  4943084 RAD cutsite not found drops (4.5%)  
 91406410 retained reads (82.7%)  
  
H9WYGADXX_1  
138847370 total sequences  
  3202776 reads contained adapter sequence (2.3%)  
   846988 barcode not found drops (0.6%)  
   267348 low quality read drops (0.2%)  
  3683206 RAD cutsite not found drops (2.7%)  
130847052 retained reads (94.2%)  
  
H9WYGADXX_2  
159871680 total sequences  
  5477072 reads contained adapter sequence (3.4%)  
  1301582 barcode not found drops (0.8%)  
   294833 low quality read drops (0.2%)  
  6884903 RAD cutsite not found drops (4.3%)  
145913290 retained reads (91.3%)  
  
HMMKLADXX_1  
147179868 total sequences  
  6490480 reads contained adapter sequence (4.4%)  
    82644 barcode not found drops (0.1%)  
   940503 low quality read drops (0.6%)  
 10192934 RAD cutsite not found drops (6.9%)  
129473307 retained reads (88.0%)  
  
HMMKLADXX_2  
149288526 total sequences  
  6056602 reads contained adapter sequence (4.1%)  
    71084 barcode not found drops (0.0%)  
  1064401 low quality read drops (0.7%)  
 12505196 RAD cutsite not found drops (8.4%)  
129591243 retained reads (86.8%)

Sum of all lanes
1,291,383,213 total sequences
67,130,654 reads contained adapter sequence (5.2%)
3,337,695 barcode not found drops (0.3%)
5,258,975 low quality read drops (0.4%)
75,209,644 RAD cutsite not found drops (5.8%)
1,140,446,245 retained reads (88.3%)

combine duplicate cleaned tags

Here we concatenate fastq files from fish with reads across multiple lanes

#!/bin/bash

# set the number of nodes
#SBATCH --nodes=1

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of output file
#SBATCH --output=

# mail alert at start, end and abortion of execution
#SBATCH --mail-type=NONE

prefixes="BP_406
BP_407
BP_408
BP_409
BP_410
BP_436
BP_437
BP_438
BP_439
BP_440
BP_466
BP_467
BP_468
BP_469
BP_470
BP_496
BP_497
BP_498
BP_499
BP_500
BP_526
BP_527
BP_528
BP_529
BP_530
BP_556
BP_557
BP_558
BP_559
BP_560
BP_586
BP_587
BP_588
BP_589
BP_590
BP_616
BP_617
BP_618
BP_619
BP_620
BP_646
BP_647
BP_648
BP_649
BP_650
BP_676
BP_677
BP_678
BP_679
BP_680
Hb_401
Hb_402
Hb_403
Hb_404
Hb_405
Hb_432
Hb_433
Hb_434
Hb_463
Hb_464
Hb_465
Hb_491
Hb_492
Hb_493
Hb_495
Hb_521
Hb_522
Hb_523
Hb_524
Hb_525
Hb_552
Hb_555
Hb_581
Hb_582
Hb_583
Hb_584
Hb_585
Hb_612
Hb_613
Hb_615
Hb_641
Hb_645
Hb_671
Hb_672
Hb_673
Hb_674
Hb_675
Mg_416
Mg_417
Mg_418
Mg_419
Mg_420
Mg_446
Mg_447
Mg_448
Mg_449
Mg_450
Mg_476
Mg_478
Mg_479
Mg_480
Mg_508
Mg_509
Mg_536
Mg_537
Mg_539
Mg_540
Mg_567
Mg_568
Mg_570
Mg_596
Mg_597
Mg_598
Mg_599
Mg_627
Mg_628
Mg_629
Mg_630
Mg_658
Mg_660
Mg_685
Mg_686
Mg_687
Mg_688
OC_421
OC_422
OC_423
OC_424
OC_425
OC_451
OC_452
OC_453
OC_454
OC_455
OC_481
OC_482
OC_483
OC_484
OC_485
OC_511
OC_512
OC_513
OC_514
OC_515
OC_541
OC_542
OC_543
OC_544
OC_545
OC_571
OC_572
OC_573
OC_574
OC_575
OC_601
OC_602
OC_603
OC_604
OC_605
OC_631
OC_632
OC_633
OC_634
OC_635
OC_661
OC_662
OC_663
OC_664
OC_665
OC_689
OC_690
RB_427
RB_428
RB_430
RB_456
RB_457
RB_458
RB_459
RB_460
RB_487
RB_488
RB_489
RB_490
RB_516
RB_519
RB_546
RB_547
RB_548
RB_549
RB_550
RB_576
RB_577
RB_579
RB_580
RB_606
RB_607
RB_608
RB_610
RB_636
RB_637
RB_638
RB_639
RB_640
RB_668
RB_669
RB_670
RB_691
RB_692
RB_693
RB_694
RB_695
SR_411
SR_412
SR_413
SR_414
SR_415
SR_441
SR_442
SR_443
SR_445
SR_471
SR_472
SR_473
SR_474
SR_475
SR_501
SR_502
SR_504
SR_531
SR_532
SR_534
SR_535
SR_562
SR_564
SR_565
SR_591
SR_594
SR_595
SR_622
SR_624
SR_625
SR_651
SR_652
SR_653
SR_655
SR_682
SR_683
SR_684
SH_2135
MG_2136
RB_2137
RB_2138
SH_2139
MG_2140
MG_2141
MG_2142
SH_2143
SH_2144
MG_2145
RB_2146
MG_2147
SH_2148
RB_2149
MG_2150
RB_2151
RB_2152
SH_2153
SH_2154
RB_2155
SH_2156
SH_2157
MG_2158
MG_2159
SH_2160
MG_2161
MG_2162
SH_2163
RB_2164
RB_2165
SH_2166
SH_2167
MG_2168
MG_2169
MG_2170
SH_2171
RB_2172
SH_2173
RB_2174
F_5001
NBH_5002
P_5003
SLC_5004
Syc_5005
H_5006
M_5007
F_5008
NBH_5009
P_5010
SLC_5011
Syc_5012
M_5013
F_5014
NBH_5015
P_5016
SLC_5017
Syc_5018
H_5019
M_5020
F_5021
NBH_5022
P_5023
SLC_5024
H_5025
M_5026
F_5027
NBH_5028
P_5029
SLC_5030
Syc_5031
H_5032
M_5033
F_5034
NBH_5035
P_5036
Syc_5037
H_5038
M_5039
F_5040
NBH_5041
P_5042
SLC_5043
Syc_5044
H_5045
M_5046
F_5047
NBH_5048
SLC_5049
Syc_5050
H_5051
M_5052
F_5053
NBH_5054
P_5055
SLC_5056
Syc_5057
H_5058
M_5059
F_5060
P_5061
SLC_5062
Syc_5063
H_5064
M_5065
F_5066
NBH_5067
P_5068
SLC_5069
Syc_5070
H_5071
M_5072
NBH_5073
P_5074
SLC_5075
Syc_5076
H_5077
M_5078
F_5079
NBH_5080
P_5081
SLC_5082
Syc_5083
H_5084
F_5085
NBH_5086
P_5087
SLC_5088
Syc_5089
H_5090
M_5091
F_5092
NBH_5093
P_5094
SLC_5095
Syc_5096
P_5097
SLC_5098
Syc_5099
H_5100
M_5101
F_5102
NBH_5103
P_5104
SLC_5105
Syc_5106
H_5107
M_5108
NBH_5109
P_5110
SLC_5111
Syc_5112
H_5113
M_5114
F_5115
NBH_5116
P_5117
SLC_5118
Syc_5119
H_5120
F_5121
NBH_5122
P_5123
SLC_5124
Syc_5125
H_5126
M_5127
F_5128
NBH_5129
P_5130
SLC_5131
Syc_5132
M_5133
F_5134
NBH_5135
P_5136
SLC_5137
Syc_5138
H_5139
M_5140
F_5141
NBH_5142
P_5143
SLC_5144
H_5145
M_5146
F_5147
NBH_5148
P_5149
SLC_5150
Syc_5151
H_5152
M_5153
F_5154
NBH_5155
P_5156
Syc_5157
H_5158
M_5159
F_5160
NBH_5161
P_5162
SLC_5163
Syc_5164
H_5165
M_5166
F_5167
NBH_5168
SLC_5169
Syc_5170
H_5171
M_5172
F_5173
NBH_5174
P_5175
SLC_5176
Syc_5177
H_5178
M_5179
F_5180
P_5181
SLC_5182
Syc_5183
H_5184
M_5185
F_5186
NBH_5187
P_5188
SLC_5189
Syc_5190
H_5191
M_5192
H_5193
M_5194
Syc_5200
H_5201
SLC_5207
Syc_5208
P_5214
SLC_5215
NBH_5221
P_5222
F_5228
NBH_5229
M_5230
M_5236
F_5237
NBH_5238
H_5244
M_5245
F_5246
" #list of duplicated sample_names - get from R

for i in $prefixes; do cat "$i"* > "$i".merged.txt ; done

for i in $prefixes; do rm "$i".fq.gz; rm "$i".1.fq.gz; done

Read Mapping

Summary

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00

#SBATCH --cpus-per-task=38

# set partition/queue to use
#SBATCH --partition=week-long-cpu

# set name of job
#SBATCH --job-name=bwa_default

# set name of output file
#SBATCH --output=bwadefault.out

# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu


source /opt/samtools-1.6/source_me

for i in `find ../cleaned_tags/ -name '*.gz'`
do
/opt/bio-bwa/bwa mem -t 38 ../genome/f_het_genome_3.0.2.fa.gz ../cleaned/${i} | /opt/samtools-1.6/bin/samtools view -@ 16 -bSu - | /opt/samtools-1.6/bin/samtools sort -@ 16 - -o ./${i:0:-6}.$

done

BAM stats

What are the depths of mapped reads (from bam files)?

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

#SBATCH --cpus-per-task=1

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of job
#SBATCH --job-name=bam_depths

# set name of output file
#SBATCH --output=bam_depths

# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu

source /opt/samtools-1.6/source_me
samtools=/opt/samtools-1.6/bin/samtools

for i in *bam
do
$samtools depth $i | awk -v var="$i" '{sum += $3; n++} END {print " '${i%.bam}' ", sum / n}' >> depths.txt
done

bam_depths<-read.table("./snp_calling/QC/depths.txt", sep = " ")
bam_depths<-bam_depths[,-c(1,3)]
colnames(bam_depths)<-c("fish_id", "depth")
ggplot(data = bam_depths)+geom_density(aes(x = depth))+labs(title = "mean depth at mapped reads (bam) over all samples", x = "depth", y = "frequency")

#we could also easily split up the depth file here to look at population level variation in depth

let’s also look at overall number mapped reads per sample

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

#SBATCH --cpus-per-task=1

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of job
#SBATCH --job-name=bam_mapped_reads

# set name of output file
#SBATCH --output=bam_mapped_reads.out


source /opt/samtools-1.6/source_me
samtools=/opt/samtools-1.6/bin/samtools

for i in *bam
do
{ $samtools view -c -F 4 $i; echo ${i%.bam} ; } | tr "\n" " " >> mapped_reads.txt
done

mapped_reads<-read.table("./snp_calling/QC/mapped_reads.txt", sep = " ")
ggplot(data = mapped_reads)+geom_histogram(aes(x = log10(mapped_reads[,1])))+labs(title =  " mapped reads in bam files", x = "log10 reads")

grandis

keeping grandis in the analysis may be useful to reconstruct the ancestral state of the genome and therefore create an unfolded SFS, but how well do the grandis reads align to the fundulus reference genome, let’s compare the bam files:

colnames(mapped_reads) <- c("reads", "sample")
mapped_reads$pop<-str_extract(mapped_reads$sample, "[^_]+")
mapped_reads %>%
  group_by(pop) %>%
  summarize(mean = mean(reads), sd = sd(reads), n = n())
ggplot()+geom_density(data = mapped_reads[(mapped_reads$pop!="grandis"),],aes(x=reads))+geom_density(data = mapped_reads[(mapped_reads$pop=="grandis"),],aes(x=reads), color = "red")+labs(title = "total number mapped reads", subtitle = "red - grandis ; black - heteroclitus")

this isn’t ideal though, lets try to figure out the PROPORTION of reads mapped to the reference genome

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

#SBATCH --cpus-per-task=1

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of job
#SBATCH --job-name=bam_mapped_reads

# set name of output file
#SBATCH --output=bam_mapped_reads.out



source /opt/samtools-1.6/source_me
samtools=/opt/samtools-1.6/bin/samtools

for i in *bam
do
{ $samtools view -c -F 0 $i; echo ${i%.bam} ; } | tr "\n" " " >> total_reads.txt
done

#total_reads <- read.table("./snp_calling/QC/total_reads.txt", sep = " ")
# make column of total reads, then proportion reads, plot again

ANGSD

ANGSD Analysis outline

We conduct two sets of ANGSD analyses. The first calculates genotype likelihoods. The second estimates the site frequency spectrum.

grViz("digraph flowchart {
      # node definitions with substituted label text for placing code
      node [fontname = Helvetica, shape = rectangle]        
      tab1 [label = '@@1']
      tab2 [label = '@@2']
      tab3 [label = '@@3']
      tab4 [label = '@@4']
      tab5 [label = '@@5']
      tab6 [label = '@@6']
      tab7 [label = '@@7']
      tab8 [label = '@@8']
      tab9 [label = '@@9']

      # edge definitions with the node IDs
      tab1 -> tab3;
      tab1 -> tab2;
      tab3 -> tab4 -> tab5
      tab4 -> tab6
      tab3 -> tab7
      tab3 -> tab8
      tab9 -> tab1
      }

      [1]: 'ANGSD -GL 1 -doMaf -doMajorMinor -glf 2'
      [2]: 'Allele Frequencies'
      [3]: 'beagle format GLs'
      [4]: 'covariance matrix (PCangsd)'
      [5]: 'PCA (R)'
      [6]: 'RDA/other multivariate EAAs (R)'
      [7]: 'LD (ngsLD)'
      [8]: 'Admixture (ngsAdmix)'
      [9]: 'bam files'
      
      ")

preservea5aff0014a782893

grViz("digraph flowchart {
      # node definitions with substituted label text
      node [fontname = Helvetica, shape = rectangle]        
      tab1 [label = '@@1']
      tab2 [label = '@@2']
      tab3 [label = '@@3']
      tab4 [label = '@@4']
      tab5 [label = '@@5']
      tab6 [label = '@@6']
      tab7 [label = '@@7']
      tab8 [label = '@@8']

      # edge definitions with the node IDs
      tab1 -> tab2
      tab2 -> tab3
      tab3 -> tab4 -> tab5
      tab3 -> tab6
      tab3 -> tab7
      tab8 -> tab1
      }

      [1]: 'ANGSD -GL 1 -doSAF unfolded'
      [2]: 'SAF'
      [3]: 'SFS (realSFS)'
      [4]: 'Dadi conversion '
      [5]: 'moments/Dadi'
      [6]: 'theta, tajD'
      [7]: 'FST'
      [8]: 'bam files'
      
      ")

preservefcb942bf6cd370cc

set up

Create lists of bams for all the separate populations


#in server directory containing bams
# for each population create a list of paths to bam files
find "$(cd ..; pwd)" -iname "BP*bam" > ../BP.bams
find "$(cd ..; pwd)" -iname "CT*bam" > ../CT.bams
find "$(cd ..; pwd)" -iname "GA*bam" > ../GA.bams
find "$(cd ..; pwd)" -iname "HB*bam" > ../HB.bams
find "$(cd ..; pwd)" -iname "KC*bam" > ../KC.bams
find "$(cd ..; pwd)" -iname "M_*bam" > ../M.bams
find "$(cd ..; pwd)" -iname "MDIBL_*bam" > ../MDIBL.bams
find "$(cd ..; pwd)" -iname "ME_*bam" > ../ME.bams
find "$(cd ..; pwd)" -iname "MG_*bam" > ../MG.bams
find "$(cd ..; pwd)" -iname "NBH_*bam" > ../NBH.bams
find "$(cd ..; pwd)" -iname "F_*bam" > ../F.bams
find "$(cd ..; pwd)" -iname "H_*bam" > ../H.bams
find "$(cd ..; pwd)" -iname "P_*bam" > ../P.bams
find "$(cd ..; pwd)" -iname "SYC_*bam" > ../SYC.bams
find "$(cd ..; pwd)" -iname "NC_*bam" > ../NC.bams
find "$(cd ..; pwd)" -iname "OC_*bam" > ../OC.bams
find "$(cd ..; pwd)" -iname "RB_*bam" > ../RB.bams
find "$(cd ..; pwd)" -iname "SC_*bam" > ../SC.bams
find "$(cd ..; pwd)" -iname "SH_*bam" > ../SH.bams
find "$(cd ..; pwd)" -iname "SLC_*bam" > ../SLC.bams
find "$(cd ..; pwd)" -iname "SR_*bam" > ../SR.bams
find "$(cd ..; pwd)" -iname "WNC_*bam" > ../WNC.bams
find "$(cd ..; pwd)" -iname "grandis_*bam" > ../grandis.bams

cat *bams > ALL.bamlist

first compress the genome and create a compression index

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

#SBATCH --cpus-per-task=10

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of job
#SBATCH --job-name=bgzip_index

# set name of output file
#SBATCH --output=bgzip.out

# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu

source /opt/samtools-1.6/source_me

/home/ddayan/software/htslib-1.9/bgzip -i -@ 10 f_het_genome_3.0.2.fa

then create a samtools index

samtools faidx f_het_genome_3.0.2.fa.gz

also need to index all the bam files

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

#SBATCH --cpus-per-task=10

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of job
#SBATCH --job-name=bgzip_index

# set name of output file
#SBATCH --output=bgzip.out

# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu

source /opt/samtools-1.6/source_me
samtools=/opt/samtools-1.6/bin/samtools

for i in ./*.bam;do $samtools index $i;done

EDA

coverage eda

First lets take a look at coverage and allele freqeuncy to get an idea what the data looks like prior to filtering

rationale for options:
command options
-P 10: use 10 threads
-ref: path to reference genome
-out: output
-r : run on only one contig (to check speed)
filtering
-uniqueOnly: use only uniquley mapping reads
-remove_bads: discard bad reads
-trim 0: perform no trimming
-C exclude reads with excessive mismatches
-baq avoid SNPs called due to indels
-minmapQ exclude reads with poor mapping
-maxDepth - r)ead depth above this level are binned (set high to llok at actual depth distrubution
dos
doqsdist calulates distribution of quality scores of used reads
dodepth calculates distrubtion of read depths
docounts needed for depth calcuation

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00

#SBATCH --cpus-per-task=10

# set partition/queue to use
#SBATCH --partition=week-long-cpu

# set name of job
#SBATCH --job-name=bwa_default

# set name of output file
#SBATCH --output=bwadefault.out

# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu

ANGSD="/opt/angsd0.928/angsd/angsd"
REF="/home/ddayan/fundulus/genome/f_het_genome_3.0.2.fa.gz"

$ANGSD -P 10 -b ../meta_data/ALL.bamlist -ref $REF -out ./Results/ALL.qc \
  -uniqueOnly 1 -remove_bads 1 -trim 0 -C 50 -baq 1 -maxDepth 1000 \
    -minMapQ 20 -doQsDist 1 -doDepth 1 -doCounts 1

results of initial ANGSD run Q scores and global depth

## barplot q-scores
qs <- read.table(file = "./snp_calling/ANGSD_outputs/qc_check/ALL.qc.qs", head=T, stringsAsFactors=F)
barplot(height=qs$counts, names.arg=qs$qscore, xlab="Q-score", ylab="Counts", main = "Q-score distribution")

## global depth
dep <- as.numeric(scan(file = "./snp_calling/ANGSD_outputs/qc_check/ALL.qc.depthGlobal",what="char", quiet=T))
dep_freq<-as.data.frame(cbind(dep, seq(1:length(dep))))
colnames(dep_freq)<-c("count", "depth")
dep_freq$depth<-dep_freq$depth-1
ggplot(data = dep_freq, aes(depth, count))+geom_bar(stat = "identity")+labs(x = "Depth", title = "global depths (summed across all fish)")

ggplot(data = dep_freq, aes(depth, count))+geom_bar(stat = "identity")+labs(x = "Depth", title = "subset of plot above")+xlim(0, 100)

Depths for a handful of individuals:

## individual depth
idep<-read.table("./snp_calling/ANGSD_outputs/qc_check/ALL.qc.depthSample")
#look at read depth histogram for the first tenindividual
i1<-as.data.frame(t(idep[1,]))
i1$depth<-c(0:999, by = 1)
colnames(i1)<-c("log10count", "depth")
i1$log10count <- log10(i1$log10count)
ggplot(data = i1, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth")

i200<-as.data.frame(t(idep[1,]))
i200$depth<-c(0:999, by = 1)
colnames(i200)<-c("log10count", "depth")
i200$log10count <- log10(i200$log10count)
ggplot(data = i200, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth")

i505<-as.data.frame(t(idep[1,]))
i505$depth<-c(0:999, by = 1)
colnames(i505)<-c("log10count", "depth")
i505$log10count <- log10(i505$log10count)
ggplot(data = i200, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth", title = )

other salient results Total number of sites analyzed: 567547 Number of sites retained after filtering: 559561

thoughts here: - most sites have very limited depth, with 1480 samples, the global depth is largely less than 1000
- 59000 sites (about 10%) have more than 1000 reads, are these driven by a moderate depth in a lot of individuals or extremely high depth in just a few?
- lets run this again but only retain sites present in 80% of individuals

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00

#SBATCH --cpus-per-task=10

# set partition/queue to use
#SBATCH --partition=week-long-cpu

# set name of job
#SBATCH --job-name=bwa_default

# set name of output file
#SBATCH --output=bwadefault.out

# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu

ANGSD="/opt/angsd0.928/angsd/angsd"
REF="/home/ddayan/fundulus/genome/f_het_genome_3.0.2.fa.gz"

$ANGSD -P 10 -b ../meta_data/ALL.bamlist -ref $REF -out ./Results/ALL.qc -r KN811289.1 \
  -uniqueOnly 1 -remove_bads 1 -trim 0 -C 50 -baq 1 -maxDepth 1000 -minInd 1184 \
    -minMapQ 20 -doQsDist 1 -doDepth 1 -doCounts 1

no sites were retained with this level of filtering, lets scale up to the ful genome and see how many sites with coverage in 80% of individuals we get: 2702 sites …

this ANGSD output is very sparse, how many loci are these sites spread across? how many reads are they based on before and after filtering? what are the mapping qualities? why is the depth data presented in the least convenient imaginable format? I miss stacks and tassel…

lets try again but a little less stringent, here we retain sites with coverage in 70% of fundulus heterociltis in the dataset (1010),

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

#SBATCH --cpus-per-task=10

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of job
#SBATCH --job-name=angsd_0.7

# set name of output file
#SBATCH --output=angsd_0.7.out


ANGSD="/opt/angsd0.928/angsd/angsd"
REF="/home/ddayan/fundulus/genome/f_het_genome_3.0.2.fa.gz"

$ANGSD -P 10 -b ../meta_data/ALL.bamlist -ref $REF -out ./Results/ALL.qc -r KN811289.1 \
  -uniqueOnly 1 -remove_bads 1 -trim 0 -C 50 -baq 1 -minInd 1010 \
    -minMapQ 20

how many loci is this over? we can use gatk countloci for this

prep:


GATK="/opt/java/jdk1.8.0_121/bin/java -jar /opt/Gatk/GenomeAnalysisTK.jar"

#create dictionary for gatk
java -jar /opt/picard-tools-1.119/CreateSequenceDictionary.jar R= fundulus/genome/f_het_genome_3.0.2.fa O= fundulus/genome/f_het_genome_3.0.2.dict
#create an index for the uncompressed genome file (gatk only run on decompressed reference genomes)
samtools faidx f_het_genome_3.0.2.fa

#now we can finally try to run GATK

$GATK -T CountLoci -R f_het_genome_3.0.2.fa -I ../aligned/OC_421.merged.bam


####agh this doesn't work because read groups were not retained in the bam files...

gl EDA (gl_0.1)

Here we’re going to finally run ANGSD

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00

#SBATCH --cpus-per-task=38

# set partition/queue to use
#SBATCH --partition=week-long-cpu

# set name of job
#SBATCH --job-name=angsd_GL_0.1

# set name of output file
#SBATCH --output=angsd_GL_0.1.out


ANGSD="/opt/angsd0.928/angsd/angsd"
REF="/home/ddayan/fundulus/genome/f_het_genome_3.0.2.fa.gz"

FILTERS="-minInd 740 -uniqueOnly 1 -remove_bads 1 -trim 0 -C 50 -baq 1 -skipTriallelic 1 -SNP_pval 1e-6 -minMapQ 20 -maxDepth 14800"
DOS="-GL 1 -doMajorMinor 1 -doMaf 2 -doGlf 2 -doCounts 1 -doDepth 1" 

$ANGSD -P 38 -b ../meta_data/ALL.bamlist -ref $REF -out ./Results/ALL.qc \
  $FILTERS $DOS \
gl_0.1 results

Total number of sites analyzed: 80148766
Number of sites retained after filtering: 62571

Coverage

## global depth
dep <- as.numeric(scan(file = "./snp_calling/ANGSD_outputs/run_0.1//ALL.qc.depthGlobal",what="char", quiet=T))
dep_freq<-as.data.frame(cbind(dep, seq(1:length(dep))))
colnames(dep_freq)<-c("count", "depth")
dep_freq$depth<-dep_freq$depth-1
dep_freq$log_count<-log10(dep_freq$count)
ggplot(data = dep_freq, aes(depth, log_count))+geom_bar(stat = "identity")+labs(x = "Depth", title = "global depths (summed across all fish)")

Of the 62571 sites, 47984 (76%) have more than 14800 reads and show up as a single thin bar at 14800 in this plot

## individual depth
idep<-read.table("./snp_calling/ANGSD_outputs/run_0.1/ALL.qc.depthSample")
#look at read depth histogram for the first tenindividual
#i100<-as.data.frame(t(idep[1,]))
#i100$depth<-c(0:14799, by = 1)
#colnames(i100)<-c("log10count", "depth")
#i100$log10count <- log10(i100$log10count)
#ggplot(data = i100, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,500))

i100<-as.data.frame(t(idep[100,]))
i100$depth<-c(0:14799, by = 1)
i100$log10count <- log10(i100$`1`)
colnames(i100)<-c("count", "depth", "log10count")
#ggplot(data = i100, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,500))
ggplot(data = i100, aes(depth, count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,400))+ylim(c(0,4000))+labs(title = "individual 100")
## Warning: Removed 14400 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (geom_bar).

i1000<-as.data.frame(t(idep[1000,]))
i1000$depth<-c(0:14799, by = 1)
i1000$log10count <- log10(i1000$`1000`)
colnames(i1000)<-c("count", "depth", "log10count")
#ggplot(data = i1000, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,1000))
ggplot(data = i1000, aes(depth, count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,400))+ylim(c(0,2000))+labs(title = "individual 1000")
## Warning: Removed 14400 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (geom_bar).

i500<-as.data.frame(t(idep[500,]))
i500$depth<-c(0:14799, by = 1)
i500$log10count <- log10(i500$`500`)
colnames(i500)<-c("count", "depth", "log10count")
#ggplot(data = i500, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,500))
ggplot(data = i500, aes(depth, count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,300))+ylim(c(0,3000))+labs(title = "individual 500")
## Warning: Removed 14500 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (geom_bar).

i148<-as.data.frame(t(idep[148,]))
i148$depth<-c(0:14799, by = 1)
i148$log10count <- log10(i148$`148`)
colnames(i148)<-c("count", "depth", "log10count")
#ggplot(data = i148, aes(depth, log10count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,148))
ggplot(data = i148, aes(depth, count))+geom_bar(stat = "identity")+labs(x = "Depth")+xlim(c(0,300))+ylim(c(0,3000))+labs(title = "grandis individual")
## Warning: Removed 14500 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (geom_bar).

if we look at a handful of depth for individual fish (reads per locus per individual) we observe a nice distribution with mostly >20x coverage

PCA eda

a note on installation (code chunk below):


#didn't want to wait for HPC to install pcangsd so did it locally. the problem (as always) is that pcangsd has depnedinecies that get install to local python but slurm calls a different python. used the tool virtualenv to get around this

#create the virtual environment in the directory containing pcagnsd
virtualenv software/pcangsd/
cd software/pcangsd/
#ctivate it
source bin/activate

#install dependcies and angsd
pip install cython
pip install numpy
python setup.py build_ext --inplace
pip install -r requirements.txt

#check install
python pcangsd.py

#deactive virtualenv

#when running in batch script make sure to source this virtual env before executing the actual command
deactivate

ok now lets try a PCA

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00

#SBATCH --cpus-per-task=38

# set partition/queue to use
#SBATCH --partition=week-long-cpu

# set name of job
#SBATCH --job-name=pcangsd_0.1

# set name of output file
#SBATCH --output=pcangsd_0.1.out



source /home/ddayan/software/pcangsd/bin/activate
python /home/ddayan/software/pcangsd/pcangsd.py -beagle ../Results/ALL_gl_0.1.beagle.gz -o pcangsd.gl_0.1 -threads 38

Let’s looks at these PCA results

cov_mat<-read.table("./pcangsd/pcangsd.gl_0.1.cov", sep = " ")

filenames = scan("./snp_calling/ANGSD_outputs/ALL.bamlist", what="character")
fish_names = gsub("/home/ddayan/fundulus/aligned/", "", gsub(".merged", "", gsub(".bam", "",  filenames)))

pc1<-prcomp(cov_mat)

pc_inds<-as.data.frame(pc1$x)
row.names(pc_inds)<-fish_names
pc_inds$pop<-gsub("_\\d+", "", fish_names)

#caution this deletes al nbh pops except nbh, fix later
pc_inds <- merge(pc_inds, sites[,c("Abbreviation", "Latitude")], by.x = "pop" , by.y ="Abbreviation")

ggplot(data = pc_inds) + geom_point(aes(x = PC1, y = PC2, color = Latitude))+scale_color_gradientn(colours = rainbow(5))#+stat_ellipse(aes(x = PC1, y = PC2, group = pop),type = "norm")

plot_ly(x=pc_inds$PC1, y=pc_inds$PC2, z=pc_inds$PC3, type="scatter3d", mode="markers", color=pc_inds$pop)

preserved7105fe5b51186fa

hover mouse over points for population id - population codes are as in the table “sites” under the section “data description” at top of document

it looks like there are some QC issues to take on (or a grandis individual is more similar to f het than other grandis), but that this approach is at least working well enough to differentiate pops in PC space

Genotype Likelihoods

filtering rationale

grandis
we will use the grandis reads to create a grandis consensus sequence (dofasta in angsd) which will be used as the ancestral state for unfolded SFS estimation, but we will not estimate genotype likelihoods in grandis as this will influence SNP probabilities and it seems that including grandis dominates the structure found by PCA

HWE
no HWE filtering, filtering out paralogues by using excess coverage instead - rationale here is that HWE is likely to be skewed given the amount of population strcutre in the sample and that filtering by HWE stratified by population is not possible in the ANGSD pipeline

MAF
MAF filtering will skew the SFS and should not be applied to the methods that rely on SFS (demography inference, FST), but the inclusion of singletons called genotyped datasets can confound STRUCTURE (little effect on multivariate population structure algorithms e.g. PCA DAPC). My assumption is that the use of genotype likelihoods instead of called genotypes should reduce the impact of singletons or low frequency alleles (erroneous or true) on the inference of structure using model based approaches like STRUCTURE, and the inclusion of rare alleles will improve discrimination among recently diverged populations (or populations with rare gene flow) when using multivariate approaches

max_coverage
will set max coverage at 2.5x the modal coverage per locus, this means will have to run ANGSD again with a higher maxdepth cutoff for the coverage output in order to estimate it…

min_coverage_locus
no hard minimum coverage cutoff, instead use stringent snp_pval cutoff and minimum number of individuals (we require that a site be present in 50% of samples - 686)

minimum_coverage_ind
the mean coverage for the individuals (at the bam level) is ~730k, but a handful of individuals have very low coverege (see bam stats section), will remove these individuals BEFORE ANGSD analysis by editing bam_lists, removed individuals in the lowest 5% of the coverage distribution

snp_pval
impose stringent probablility cutoff for reatining a SNP: 1e-6

quality_and_mapping
-uniqueOnly 1 -remove_bads 1 -C 50 -baq 1 -skipTriallelic 1 -minMapQ 20

ANGSD GL Run 1.2x

version history 1.0 first run 1.1 removed outlier high coverage 1.2 start from scratch with mislabeled GA fish removed - GA_3318 clustered with grandis in PC space and also was it’s own group in PCA run on genotype likelihoods v1.1 1.3 coverage filtered output of 1.2

Make list of bam files, filtering out grandis and individuals in the lowest 5% of reads:

#get rid of individuals with low sequencing
bad_inds <- mapped_reads[mapped_reads$reads<quantile(mapped_reads$reads, 0.05),2]
#also get rid of mislabeled individual - pca_1.1 showed that "GA_3318" is most likely a grandis
bad_inds<-append(bad_inds, mapped_reads[140,2])
bad_inds <- mapped_reads[bad_inds,2]
#write.table(bad_inds, file = "./snp_calling/QC/bad_inds.txt", quote = FALSE, row.names = FALSE)

grep -v -f ./snp_calling/QC/bad_inds.txt ./snp_calling/ANGSD_outputs/ALL.bamlist | grep -v 'grandis' > ./snp_calling/QC/good_no_grandis_2.bamlist

This removed all grandis, 72 individuals with too few reads and 1 individual labeled as GA but actually grandis, leaving a final dataset of 1371 individuals.

Make genotype likelihood estimation v1.2

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00

#SBATCH --cpus-per-task=38

# set partition/queue to use
#SBATCH --partition=week-long-cpu

# set name of job
#SBATCH --job-name=angsd_GL_1.0

# set name of output file
#SBATCH --output=angsd_GL_1.0.out


ANGSD="/opt/angsd0.928/angsd/angsd"
REF="/home/ddayan/fundulus/genome/f_het_genome_3.0.2.fa.gz"

FILTERS="-minInd 686 -uniqueOnly 1 -remove_bads 1 -trim 0 -C 50 -baq 1 -skipTriallelic 1 -SNP_pval 1e-6 -minMapQ 20"
DOS="-GL 1 -doMajorMinor 1 -doMaf 2 -doGlf 2 -doCounts 1 -dumpCounts 2" 

$ANGSD -P 38 -b ../meta_data/good_no_grandis_2.bamlist -ref $REF -out ./Results/gl_1.2 \
  $FILTERS $DOS \

Run 1.0: Total number of sites analyzed: 78663923 Number of sites retained after filtering: 77717

Run 1.2 Total number of sites analyzed: 78419827 Number of sites retained after filtering: 75430

Coverage

Extract and plot coverage data

#check the format before running this
depths_loci<-read.table("./snp_calling/QC/gl_1.2.pos", header = T, sep = "\t")

#mode calculation
getMode <- function(x) {
keys <- unique(x)
keys[which.max(tabulate(match(x, keys)))]
}

modal_depth <- getMode(depths_loci$totDepth)
ggplot(data = depths_loci) + geom_histogram(aes(x = totDepth))+xlim(c(0,75000)) + geom_vline(aes(xintercept = modal_depth), color = "red") + geom_vline(aes(xintercept = 2.5*modal_depth), color = "blue")

The modal depth per locus is modal_depth. All loci with 2.5x greater read depth are then filtered out.

#note, remove comments if need to write the outputs to file

#find loci with >2.5x the mode and filter then write out to loci list which will be used to filter the GL outputs (filtered datasets are named 1.1)
sites_to_keep<-subset(depths_loci, totDepth < 2.5*modal_depth)
#write.table(sites_to_keep[,c(1,2)], file = "./snp_calling/QC/sites_to_keep.txt", row.names = FALSE, quote = FALSE, col.names = FALSE)

#get bad sites and store in beagle call format
sites_to_remove<-subset(depths_loci, totDepth > 2.5*modal_depth)
sites_to_remove$beagle_id<-paste(sites_to_remove$chr, "_", sites_to_remove$pos, sep = "")
#write.table(sites_to_remove[,c(4)], file = "./snp_calling/QC/sites_to_remove.txt", row.names = FALSE, quote = FALSE, col.names = FALSE)

GL_1.2 had nrow(depths_loci) SNPs, after fitlering out SNPs with greater than 2.5x the mode coverage (modal_depthx = mode) there are nrow(sites_to_keep) SNPs

#filtering the datasets

zcat gl_1.2.beagle.gz > gl_1.2.beagle 
awk 'FNR==NR { a[$1]; next } !($1 in a){print $0}'  sites_to_remove.txt gl_1.2.beagle  > gl_1.3.beagle
bgzip gl_1.3.beagle

Population Structure

PCA

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00

#SBATCH --cpus-per-task=38

# set partition/queue to use
#SBATCH --partition=week-long-cpu

# set name of job
#SBATCH --job-name=pcangsd_1.3

# set name of output file
#SBATCH --output=pcangsd_1.3.out



source /home/ddayan/software/pcangsd/bin/activate
python /home/ddayan/software/pcangsd/pcangsd.py -beagle ../Results/gl_1.3.beagle.gz -o pcangsd.gl_1.1 -threads 38 -minMaf 0.0

PCANGSD run info:

Read 1371 samples and 73375 sites

Estimating population allele frequencies EM (MAF) converged at iteration: 12

Estimating covariance matrix Using 5 principal components (MAP test)

cov_mat<-read.table("./pcangsd/pcangsd.gl_1.3.cov", sep = " ")

filenames <- mapped_reads[! mapped_reads$sample %in% bad_inds,]
filenames <- filenames[!grepl("grandis", filenames$sample),]
filenames <- filenames[-c(1372, 1373),]


fish_names =  gsub(".merged", "", filenames$sample)

pc1<-prcomp(cov_mat)

save(pc1, file = "./pcangsd/pca.R")

pc_inds<-as.data.frame(pc1$x)
row.names(pc_inds)<-fish_names
pc_inds$pop<-gsub("_\\d+", "", fish_names)

# write.table(pc_inds, file = "./pcangsd/pca.txt", quote = FALSE) #commented to avoid overwriting, run if doing first run through analysis


#caution this deletes al nbh pops except nbh, fix later
pc_inds <- merge(pc_inds, sites[,c("Abbreviation", "Latitude")], by.x = "pop" , by.y ="Abbreviation")

ggplot(data = pc_inds) + geom_point(aes(x = PC1, y = PC2, color = Latitude))+scale_color_gradientn(colours = rainbow(5))#+stat_ellipse(aes(x = PC1, y = PC2, group = pop),type = "norm")

plot_ly(x=pc_inds$PC1, y=pc_inds$PC2, z=pc_inds$PC3, type="scatter3d", mode="markers", color=pc_inds$pop)

preserve506392caba688e5e

There is 1 outlier fish (“GA_3318”). Given that this fish was caught in Georgia, it is likely a misidentified grandis is needs to be excluded from the analysis. New GL file and all resulting analysis labeled 1.2

Admixture

NGSadmix is used to estimate ancestry for K = 1-6 with 10 replicate runs and best k determined by evanno

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00

#SBATCH --cpus-per-task=10

# set partition/queue to use
#SBATCH --partition=week-long-cpu

# set name of job
#SBATCH --job-name=ngsadmix_1.3_k9

# set name of output file
#SBATCH --output=ngsadmix_1.3_k9


NGSadmix="/opt/angsd0.928/angsd/misc/NGSadmix"
BEAGLE="/home/ddayan/fundulus/ANGSD/Results/gl_1.3.beagle.gz"

Out_Prefix="1.3_k9"
K="9"

$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.1 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.2 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.3 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.4 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.5 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.6 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.7 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.8 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.9 -P 10 -minMaf 0
$NGSadmix -likes $BEAGLE -K $K -o $Out_Prefix.10 -P 10 -minMaf 0

Plot admixture portions for a single run of k = 2

#make popfile
pops<- as.data.frame(cbind(fish_names, gsub("_\\d+", "", fish_names)))
colnames(pops) <- c("sample", "pop")

#plot a admixture plot
q<-read.table("./ngs_admix/sandbox/outFileName.qopt")
    #set order of poulations by latitude

q$sample <- pops$sample
q$pop <- pops$pop

pops$pop <- factor( as.character(pops$pop), levels=c("GA", "SC", "WNC", "NC", "KC", "SH", "RB", "OC", "MG", "Mg", "CT", "SR", "BP", "HB", "Hb", "H", "F", "M", "Syc", "NBH", "P", "SLC", "ME", "MDIBL") )
ord <- order(pops$pop)

q <- q[ord,]

plot_data <- q %>% 
  mutate(id = sample) %>% 
  gather('cluster', 'prob', V1:V2) %>% 
  group_by(id) %>% 
  mutate(likely_assignment = cluster[which.max(prob)],
         assingment_prob = max(prob)) %>% 
  arrange(likely_assignment, desc(assingment_prob)) %>% 
  ungroup() %>% 
  dplyr::arrange(factor( as.character(pop), levels=c("GA", "SC", "WNC", "NC", "KC", "SH", "RB", "OC", "MG", "Mg", "CT", "SR", "BP", "HB", "Hb", "H", "F", "M", "Syc", "NBH", "P", "SLC", "ME", "MDIBL")  ))

#reorder levels
plot_data$pop <- factor( as.character(plot_data$pop), levels=c("GA", "SC", "WNC", "NC", "KC", "SH", "RB", "OC", "MG", "Mg", "CT", "SR", "BP", "HB", "Hb", "H", "F", "M", "Syc", "NBH", "P", "SLC", "ME", "MDIBL")  )



#ggplot(plot_data, aes(id, prob, fill = cluster)) +
#  geom_col() +
#  theme_classic()

ggplot(plot_data, aes(id, prob, fill = cluster)) +
  geom_col() +
  facet_grid(~pop, scales = 'free', space = 'free')+ theme(panel.spacing=unit(0.01, "lines"), axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank(), legend.position = "none")

format data for clumpak


(for log in `ls *.log`; do grep -Po 'like=\K[^ ]+' $log; done) > logfile

logs <- as.data.frame(read.table("./ngs_admix/logfile_likelihoods.txt"))
logs$K <- c(rep("1", 10), rep("2", 
    10), rep("3", 10), rep("4", 10), rep("5", 10), rep("6", 10))
colnames(logs)[1] <- "likelihood"
write.table(logs[, c(2, 1)], "/ngs_admix/logfile_likelihood_formatted.txt", row.names = F, 
    col.names = F, quote = F)

delta K method suggest best K is 2, but prob(k) is maximixed at k = 6, therefore will run for more K

SFS

SFS Analysis Outline
* Make consensus sequence from grandis reads to polarize the SAF - doFasta -2
* estimate SAFs for each population and metapopulation -doSAF * estimate per population SFS
* estimate 2dSFS for each pair of populations
* 1d and 2dSFS used in downstream analyses: per-site summary statistics, global Fst, and dadi

questions: is it possible to estimate SAF without re-estimating GLs (we already have them) what filtering is neccesarry for the population and 2d sfs, when there are so few indiviudals a maf of 1% becomes unmeaningful, should we filter to a higher MAF or at least, increase min number of individuals for a site to be retained (i.e. saf estimate is likely to be very very bad if there are sites with 2 individuals) *how should we pool sampling locales into metapopulations (i.e. it’s clear that all the NBH populations should be combined into a single pop, but what about all NJ populations or all ME populations) - what spatial or genetic scale should be the cutoff for pooling populations

1d SFS estimation

ancestral sequence

Polarize SAF by an outgroup (f grandis) use doFasta -2 with -doCounts 1 and the same filters as $FILTERS from the GL run


#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00

#SBATCH --cpus-per-task=38

# set partition/queue to use
#SBATCH --partition=week-long-cpu

# set name of job
#SBATCH --job-name=consensus

# set name of output file
#SBATCH --output=consensus.out


ANGSD="/opt/angsd0.928/angsd/angsd"
REF="/home/ddayan/fundulus/genome/f_het_genome_3.0.2.fa.gz"

FILTERS="-uniqueOnly 1 -remove_bads 1 -trim 0 -C 50 -baq 1 -minMapQ 20"
DOS="-doFasta 2 -doCounts 1 -dumpCounts 2" 

$ANGSD -P 38 -b ../meta_data/grandis.bams -ref $REF -out ./Results/grandis_consensus.fa \
  $FILTERS $DOS \

SAF 1 big population

2d SFS

global FST

Per site summary stats

Landscape Genetics

require(vegan)
require(ade4)
require(adespatial)
require(SoDA)
require(marmap)

spatial and environmental data

distances
rather than using distance matrices based on latitude and longitude, we make a more biologically relevant model where spatial distances among populations is determined by connections via shallow water <5 meters) this is accomplished via the marmap package in R

atlantic <- getNOAA.bathy(lon1 = -83, lon2 = -65, lat1 = 28, lat2 = 48, resolution = 1)
blues <- c("lightsteelblue4", "lightsteelblue3", "lightsteelblue2", "lightsteelblue1")
greys <- c(grey(0.6), grey(0.93), grey(0.99))
plot(atlantic, image = TRUE, land = TRUE, lwd = 0.03, bpal = list(c(0, max(atlantic), greys),c(min(atlantic), 0, blues)))
#sites<-read.table("~/Science/fundulus_clinal/cline/environmental/lat_lon2.txt", header=T)
sites<-read.table("./environmental/bathy_sites.txt", header = T, sep = "\t")
#remove braytonpt and oc
#sites3<-sites2[c(-7,-11),]
#row.names(sites)<-sites[,1]
#sites<-sites[,-1]
#sites2<-sites
#sites2[,1]<-sites[,2]
#sites2[,2]<-sites[,1]
trans<-trans.mat(atlantic, min.depth=1 , max.depth=-15)
path <- lc.dist(trans, sites[,c(3,2)], res = "path")
##          lon      lat depth
## 1  -71.18604 41.71250    -5
## 2  -72.52448 41.27902     0
## 3  -81.35452 31.44932    -2
## 4  -71.03655 41.49594    -3
## 5  -76.42462 37.29845     3
## 6  -70.81464 41.65875    -5
## 7  -68.31693 44.44181    41
## 8  -69.72242 43.94622     0
## 9  -74.07045 40.05046    -1
## 10 -70.90760 41.62486    -6
## 11 -75.65578 35.91287     0
## 12 -74.18475 39.80866     1
## 13 -74.32448 39.50878    -3
## 14 -79.91624 32.75873    -5
## 15 -74.75746 39.06786    -3
## 16 -70.97979 41.53583    -1
## 17 -71.52557 41.38235    -6
## 18 -77.92477 34.02835     0
## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |                                                                 |   1%
  |                                                                       
  |=                                                                |   1%
  |                                                                       
  |=                                                                |   2%
  |                                                                       
  |==                                                               |   3%
  |                                                                       
  |===                                                              |   4%
  |                                                                       
  |===                                                              |   5%
  |                                                                       
  |====                                                             |   6%
  |                                                                       
  |====                                                             |   7%
  |                                                                       
  |=====                                                            |   7%
  |                                                                       
  |=====                                                            |   8%
  |                                                                       
  |======                                                           |   8%
  |                                                                       
  |======                                                           |   9%
  |                                                                       
  |======                                                           |  10%
  |                                                                       
  |=======                                                          |  10%
  |                                                                       
  |=======                                                          |  11%
  |                                                                       
  |========                                                         |  12%
  |                                                                       
  |========                                                         |  13%
  |                                                                       
  |=========                                                        |  14%
  |                                                                       
  |==========                                                       |  15%
  |                                                                       
  |==========                                                       |  16%
  |                                                                       
  |===========                                                      |  16%
  |                                                                       
  |===========                                                      |  17%
  |                                                                       
  |===========                                                      |  18%
  |                                                                       
  |============                                                     |  18%
  |                                                                       
  |============                                                     |  19%
  |                                                                       
  |=============                                                    |  20%
  |                                                                       
  |==============                                                   |  21%
  |                                                                       
  |==============                                                   |  22%
  |                                                                       
  |===============                                                  |  23%
  |                                                                       
  |===============                                                  |  24%
  |                                                                       
  |================                                                 |  24%
  |                                                                       
  |================                                                 |  25%
  |                                                                       
  |=================                                                |  25%
  |                                                                       
  |=================                                                |  26%
  |                                                                       
  |=================                                                |  27%
  |                                                                       
  |==================                                               |  27%
  |                                                                       
  |==================                                               |  28%
  |                                                                       
  |===================                                              |  29%
  |                                                                       
  |====================                                             |  30%
  |                                                                       
  |====================                                             |  31%
  |                                                                       
  |=====================                                            |  32%
  |                                                                       
  |=====================                                            |  33%
  |                                                                       
  |======================                                           |  33%
  |                                                                       
  |======================                                           |  34%
  |                                                                       
  |=======================                                          |  35%
  |                                                                       
  |=======================                                          |  36%
  |                                                                       
  |========================                                         |  37%
  |                                                                       
  |=========================                                        |  38%
  |                                                                       
  |=========================                                        |  39%
  |                                                                       
  |==========================                                       |  40%
  |                                                                       
  |==========================                                       |  41%
  |                                                                       
  |===========================                                      |  41%
  |                                                                       
  |===========================                                      |  42%
  |                                                                       
  |============================                                     |  42%
  |                                                                       
  |============================                                     |  43%
  |                                                                       
  |============================                                     |  44%
  |                                                                       
  |=============================                                    |  44%
  |                                                                       
  |=============================                                    |  45%
  |                                                                       
  |==============================                                   |  46%
  |                                                                       
  |===============================                                  |  47%
  |                                                                       
  |===============================                                  |  48%
  |                                                                       
  |================================                                 |  49%
  |                                                                       
  |================================                                 |  50%
  |                                                                       
  |=================================                                |  50%
  |                                                                       
  |=================================                                |  51%
  |                                                                       
  |==================================                               |  52%
  |                                                                       
  |==================================                               |  53%
  |                                                                       
  |===================================                              |  54%
  |                                                                       
  |====================================                             |  55%
  |                                                                       
  |====================================                             |  56%
  |                                                                       
  |=====================================                            |  56%
  |                                                                       
  |=====================================                            |  57%
  |                                                                       
  |=====================================                            |  58%
  |                                                                       
  |======================================                           |  58%
  |                                                                       
  |======================================                           |  59%
  |                                                                       
  |=======================================                          |  59%
  |                                                                       
  |=======================================                          |  60%
  |                                                                       
  |========================================                         |  61%
  |                                                                       
  |========================================                         |  62%
  |                                                                       
  |=========================================                        |  63%
  |                                                                       
  |==========================================                       |  64%
  |                                                                       
  |==========================================                       |  65%
  |                                                                       
  |===========================================                      |  66%
  |                                                                       
  |===========================================                      |  67%
  |                                                                       
  |============================================                     |  67%
  |                                                                       
  |============================================                     |  68%
  |                                                                       
  |=============================================                    |  69%
  |                                                                       
  |=============================================                    |  70%
  |                                                                       
  |==============================================                   |  71%
  |                                                                       
  |===============================================                  |  72%
  |                                                                       
  |===============================================                  |  73%
  |                                                                       
  |================================================                 |  73%
  |                                                                       
  |================================================                 |  74%
  |                                                                       
  |================================================                 |  75%
  |                                                                       
  |=================================================                |  75%
  |                                                                       
  |=================================================                |  76%
  |                                                                       
  |==================================================               |  76%
  |                                                                       
  |==================================================               |  77%
  |                                                                       
  |===================================================              |  78%
  |                                                                       
  |===================================================              |  79%
  |                                                                       
  |====================================================             |  80%
  |                                                                       
  |=====================================================            |  81%
  |                                                                       
  |=====================================================            |  82%
  |                                                                       
  |======================================================           |  82%
  |                                                                       
  |======================================================           |  83%
  |                                                                       
  |======================================================           |  84%
  |                                                                       
  |=======================================================          |  84%
  |                                                                       
  |=======================================================          |  85%
  |                                                                       
  |========================================================         |  86%
  |                                                                       
  |=========================================================        |  87%
  |                                                                       
  |=========================================================        |  88%
  |                                                                       
  |==========================================================       |  89%
  |                                                                       
  |==========================================================       |  90%
  |                                                                       
  |===========================================================      |  90%
  |                                                                       
  |===========================================================      |  91%
  |                                                                       
  |===========================================================      |  92%
  |                                                                       
  |============================================================     |  92%
  |                                                                       
  |============================================================     |  93%
  |                                                                       
  |=============================================================    |  93%
  |                                                                       
  |=============================================================    |  94%
  |                                                                       
  |==============================================================   |  95%
  |                                                                       
  |==============================================================   |  96%
  |                                                                       
  |===============================================================  |  97%
  |                                                                       
  |================================================================ |  98%
  |                                                                       
  |================================================================ |  99%
  |                                                                       
  |=================================================================|  99%
  |                                                                       
  |=================================================================| 100%
#distance_matrix <- lc.dist(trans, sites[,c(3,2)], res = "dist")
points(sites[,c(3,2)], pch = 21, col = "blue", bg = col2alpha("blue", .9),cex = 1.2)
lapply(path, lines, col = "orange", lwd = 2, lty = 1)->dummy

#get.depth(atlantic, sites[,c(3,2)] , locator = FALSE)

#distance_matrix
#prep site data for dbmem
  #merge all NBH inds to a single pop and convert to all caps
pops <- pops %>%  
  mutate(pop = as.character(pop)) %>%
  mutate(pop = replace(pop, pop == "F", "NBH")) %>%
  mutate(pop = replace(pop, pop == "P" , "NBH")) %>%
  mutate(pop = replace(pop, pop ==  "Syc" , "NBH")) %>%
  mutate(pop = replace(pop, pop ==  "H", "NBH")) %>%
  mutate(pop = replace(pop, pop == "Mg", "MG")) %>%
  mutate(pop = replace(pop, pop == "Hb", "HB"))


ind_sites <- pops %>%
  left_join(., sites, by = "pop")
## Warning: Column `pop` joining character vector and factor, coercing into
## character vector
distance_matrix <- lc.dist(trans, ind_sites[,c(4,3)], res = "dist")
##            lon      lat depth
## 1    -71.18604 41.71250    -5
## 2    -71.18604 41.71250    -5
## 3    -71.18604 41.71250    -5
## 4    -71.18604 41.71250    -5
## 5    -71.18604 41.71250    -5
## 6    -71.18604 41.71250    -5
## 7    -71.18604 41.71250    -5
## 8    -71.18604 41.71250    -5
## 9    -71.18604 41.71250    -5
## 10   -71.18604 41.71250    -5
## 11   -71.18604 41.71250    -5
## 12   -71.18604 41.71250    -5
## 13   -71.18604 41.71250    -5
## 14   -71.18604 41.71250    -5
## 15   -71.18604 41.71250    -5
## 16   -71.18604 41.71250    -5
## 17   -71.18604 41.71250    -5
## 18   -71.18604 41.71250    -5
## 19   -71.18604 41.71250    -5
## 20   -71.18604 41.71250    -5
## 21   -71.18604 41.71250    -5
## 22   -71.18604 41.71250    -5
## 23   -71.18604 41.71250    -5
## 24   -71.18604 41.71250    -5
## 25   -71.18604 41.71250    -5
## 26   -71.18604 41.71250    -5
## 27   -71.18604 41.71250    -5
## 28   -71.18604 41.71250    -5
## 29   -71.18604 41.71250    -5
## 30   -71.18604 41.71250    -5
## 31   -71.18604 41.71250    -5
## 32   -71.18604 41.71250    -5
## 33   -71.18604 41.71250    -5
## 34   -71.18604 41.71250    -5
## 35   -71.18604 41.71250    -5
## 36   -71.18604 41.71250    -5
## 37   -71.18604 41.71250    -5
## 38   -71.18604 41.71250    -5
## 39   -71.18604 41.71250    -5
## 40   -71.18604 41.71250    -5
## 41   -71.18604 41.71250    -5
## 42   -71.18604 41.71250    -5
## 43   -71.18604 41.71250    -5
## 44   -71.18604 41.71250    -5
## 45   -71.18604 41.71250    -5
## 46   -71.18604 41.71250    -5
## 47   -71.18604 41.71250    -5
## 48   -71.18604 41.71250    -5
## 49   -71.18604 41.71250    -5
## 50   -71.18604 41.71250    -5
## 51   -72.52448 41.27902     0
## 52   -72.52448 41.27902     0
## 53   -72.52448 41.27902     0
## 54   -72.52448 41.27902     0
## 55   -72.52448 41.27902     0
## 56   -72.52448 41.27902     0
## 57   -72.52448 41.27902     0
## 58   -72.52448 41.27902     0
## 59   -72.52448 41.27902     0
## 60   -72.52448 41.27902     0
## 61   -72.52448 41.27902     0
## 62   -72.52448 41.27902     0
## 63   -72.52448 41.27902     0
## 64   -72.52448 41.27902     0
## 65   -72.52448 41.27902     0
## 66   -72.52448 41.27902     0
## 67   -72.52448 41.27902     0
## 68   -72.52448 41.27902     0
## 69   -72.52448 41.27902     0
## 70   -72.52448 41.27902     0
## 71   -70.90760 41.62486    -6
## 72   -70.90760 41.62486    -6
## 73   -70.90760 41.62486    -6
## 74   -70.90760 41.62486    -6
## 75   -70.90760 41.62486    -6
## 76   -70.90760 41.62486    -6
## 77   -70.90760 41.62486    -6
## 78   -70.90760 41.62486    -6
## 79   -70.90760 41.62486    -6
## 80   -70.90760 41.62486    -6
## 81   -70.90760 41.62486    -6
## 82   -70.90760 41.62486    -6
## 83   -70.90760 41.62486    -6
## 84   -70.90760 41.62486    -6
## 85   -70.90760 41.62486    -6
## 86   -70.90760 41.62486    -6
## 87   -70.90760 41.62486    -6
## 88   -70.90760 41.62486    -6
## 89   -70.90760 41.62486    -6
## 90   -70.90760 41.62486    -6
## 91   -70.90760 41.62486    -6
## 92   -70.90760 41.62486    -6
## 93   -70.90760 41.62486    -6
## 94   -70.90760 41.62486    -6
## 95   -70.90760 41.62486    -6
## 96   -70.90760 41.62486    -6
## 97   -70.90760 41.62486    -6
## 98   -70.90760 41.62486    -6
## 99   -70.90760 41.62486    -6
## 100  -70.90760 41.62486    -6
## 101  -81.35452 31.44932    -2
## 102  -81.35452 31.44932    -2
## 103  -81.35452 31.44932    -2
## 104  -81.35452 31.44932    -2
## 105  -81.35452 31.44932    -2
## 106  -81.35452 31.44932    -2
## 107  -81.35452 31.44932    -2
## 108  -81.35452 31.44932    -2
## 109  -81.35452 31.44932    -2
## 110  -81.35452 31.44932    -2
## 111  -81.35452 31.44932    -2
## 112  -81.35452 31.44932    -2
## 113  -81.35452 31.44932    -2
## 114  -81.35452 31.44932    -2
## 115  -81.35452 31.44932    -2
## 116  -81.35452 31.44932    -2
## 117  -81.35452 31.44932    -2
## 118  -81.35452 31.44932    -2
## 119  -81.35452 31.44932    -2
## 120  -81.35452 31.44932    -2
## 121  -81.35452 31.44932    -2
## 122  -81.35452 31.44932    -2
## 123  -81.35452 31.44932    -2
## 124  -81.35452 31.44932    -2
## 125  -81.35452 31.44932    -2
## 126  -81.35452 31.44932    -2
## 127  -81.35452 31.44932    -2
## 128  -81.35452 31.44932    -2
## 129  -81.35452 31.44932    -2
## 130  -81.35452 31.44932    -2
## 131  -81.35452 31.44932    -2
## 132  -81.35452 31.44932    -2
## 133  -81.35452 31.44932    -2
## 134  -81.35452 31.44932    -2
## 135  -81.35452 31.44932    -2
## 136  -81.35452 31.44932    -2
## 137  -81.35452 31.44932    -2
## 138  -81.35452 31.44932    -2
## 139  -70.90760 41.62486    -6
## 140  -70.90760 41.62486    -6
## 141  -70.90760 41.62486    -6
## 142  -70.90760 41.62486    -6
## 143  -70.90760 41.62486    -6
## 144  -70.90760 41.62486    -6
## 145  -70.90760 41.62486    -6
## 146  -70.90760 41.62486    -6
## 147  -70.90760 41.62486    -6
## 148  -70.90760 41.62486    -6
## 149  -70.90760 41.62486    -6
## 150  -70.90760 41.62486    -6
## 151  -70.90760 41.62486    -6
## 152  -70.90760 41.62486    -6
## 153  -70.90760 41.62486    -6
## 154  -70.90760 41.62486    -6
## 155  -70.90760 41.62486    -6
## 156  -70.90760 41.62486    -6
## 157  -70.90760 41.62486    -6
## 158  -70.90760 41.62486    -6
## 159  -70.90760 41.62486    -6
## 160  -70.90760 41.62486    -6
## 161  -70.90760 41.62486    -6
## 162  -70.90760 41.62486    -6
## 163  -70.90760 41.62486    -6
## 164  -70.90760 41.62486    -6
## 165  -70.90760 41.62486    -6
## 166  -70.90760 41.62486    -6
## 167  -70.90760 41.62486    -6
## 168  -70.90760 41.62486    -6
## 169  -71.03655 41.49594    -3
## 170  -71.03655 41.49594    -3
## 171  -71.03655 41.49594    -3
## 172  -71.03655 41.49594    -3
## 173  -71.03655 41.49594    -3
## 174  -71.03655 41.49594    -3
## 175  -71.03655 41.49594    -3
## 176  -71.03655 41.49594    -3
## 177  -71.03655 41.49594    -3
## 178  -71.03655 41.49594    -3
## 179  -71.03655 41.49594    -3
## 180  -71.03655 41.49594    -3
## 181  -71.03655 41.49594    -3
## 182  -71.03655 41.49594    -3
## 183  -71.03655 41.49594    -3
## 184  -71.03655 41.49594    -3
## 185  -71.03655 41.49594    -3
## 186  -71.03655 41.49594    -3
## 187  -71.03655 41.49594    -3
## 188  -71.03655 41.49594    -3
## 189  -71.03655 41.49594    -3
## 190  -71.03655 41.49594    -3
## 191  -71.03655 41.49594    -3
## 192  -71.03655 41.49594    -3
## 193  -71.03655 41.49594    -3
## 194  -71.03655 41.49594    -3
## 195  -71.03655 41.49594    -3
## 196  -71.03655 41.49594    -3
## 197  -71.03655 41.49594    -3
## 198  -71.03655 41.49594    -3
## 199  -71.03655 41.49594    -3
## 200  -71.03655 41.49594    -3
## 201  -71.03655 41.49594    -3
## 202  -71.03655 41.49594    -3
## 203  -71.03655 41.49594    -3
## 204  -71.03655 41.49594    -3
## 205  -71.03655 41.49594    -3
## 206  -71.03655 41.49594    -3
## 207  -71.03655 41.49594    -3
## 208  -71.03655 41.49594    -3
## 209  -71.03655 41.49594    -3
## 210  -71.03655 41.49594    -3
## 211  -71.03655 41.49594    -3
## 212  -71.03655 41.49594    -3
## 213  -71.03655 41.49594    -3
## 214  -71.03655 41.49594    -3
## 215  -71.03655 41.49594    -3
## 216  -71.03655 41.49594    -3
## 217  -76.42462 37.29845     3
## 218  -76.42462 37.29845     3
## 219  -76.42462 37.29845     3
## 220  -76.42462 37.29845     3
## 221  -76.42462 37.29845     3
## 222  -76.42462 37.29845     3
## 223  -76.42462 37.29845     3
## 224  -76.42462 37.29845     3
## 225  -76.42462 37.29845     3
## 226  -76.42462 37.29845     3
## 227  -76.42462 37.29845     3
## 228  -76.42462 37.29845     3
## 229  -76.42462 37.29845     3
## 230  -76.42462 37.29845     3
## 231  -76.42462 37.29845     3
## 232  -76.42462 37.29845     3
## 233  -76.42462 37.29845     3
## 234  -76.42462 37.29845     3
## 235  -76.42462 37.29845     3
## 236  -76.42462 37.29845     3
## 237  -76.42462 37.29845     3
## 238  -76.42462 37.29845     3
## 239  -76.42462 37.29845     3
## 240  -76.42462 37.29845     3
## 241  -76.42462 37.29845     3
## 242  -70.81464 41.65875    -5
## 243  -70.81464 41.65875    -5
## 244  -70.81464 41.65875    -5
## 245  -70.81464 41.65875    -5
## 246  -70.81464 41.65875    -5
## 247  -70.81464 41.65875    -5
## 248  -70.81464 41.65875    -5
## 249  -70.81464 41.65875    -5
## 250  -70.81464 41.65875    -5
## 251  -70.81464 41.65875    -5
## 252  -70.81464 41.65875    -5
## 253  -70.81464 41.65875    -5
## 254  -70.81464 41.65875    -5
## 255  -70.81464 41.65875    -5
## 256  -70.81464 41.65875    -5
## 257  -70.81464 41.65875    -5
## 258  -70.81464 41.65875    -5
## 259  -70.81464 41.65875    -5
## 260  -70.81464 41.65875    -5
## 261  -70.81464 41.65875    -5
## 262  -70.81464 41.65875    -5
## 263  -70.81464 41.65875    -5
## 264  -70.81464 41.65875    -5
## 265  -70.81464 41.65875    -5
## 266  -70.81464 41.65875    -5
## 267  -70.81464 41.65875    -5
## 268  -70.81464 41.65875    -5
## 269  -70.81464 41.65875    -5
## 270  -70.81464 41.65875    -5
## 271  -70.81464 41.65875    -5
## 272  -70.81464 41.65875    -5
## 273  -68.31693 44.44181    41
## 274  -68.31693 44.44181    41
## 275  -68.31693 44.44181    41
## 276  -68.31693 44.44181    41
## 277  -68.31693 44.44181    41
## 278  -68.31693 44.44181    41
## 279  -68.31693 44.44181    41
## 280  -68.31693 44.44181    41
## 281  -68.31693 44.44181    41
## 282  -68.31693 44.44181    41
## 283  -68.31693 44.44181    41
## 284  -68.31693 44.44181    41
## 285  -68.31693 44.44181    41
## 286  -68.31693 44.44181    41
## 287  -68.31693 44.44181    41
## 288  -68.31693 44.44181    41
## 289  -68.31693 44.44181    41
## 290  -68.31693 44.44181    41
## 291  -68.31693 44.44181    41
## 292  -68.31693 44.44181    41
## 293  -68.31693 44.44181    41
## 294  -68.31693 44.44181    41
## 295  -68.31693 44.44181    41
## 296  -69.72242 43.94622     0
## 297  -69.72242 43.94622     0
## 298  -69.72242 43.94622     0
## 299  -69.72242 43.94622     0
## 300  -69.72242 43.94622     0
## 301  -69.72242 43.94622     0
## 302  -69.72242 43.94622     0
## 303  -69.72242 43.94622     0
## 304  -69.72242 43.94622     0
## 305  -69.72242 43.94622     0
## 306  -69.72242 43.94622     0
## 307  -69.72242 43.94622     0
## 308  -69.72242 43.94622     0
## 309  -69.72242 43.94622     0
## 310  -69.72242 43.94622     0
## 311  -69.72242 43.94622     0
## 312  -69.72242 43.94622     0
## 313  -69.72242 43.94622     0
## 314  -69.72242 43.94622     0
## 315  -69.72242 43.94622     0
## 316  -69.72242 43.94622     0
## 317  -69.72242 43.94622     0
## 318  -69.72242 43.94622     0
## 319  -69.72242 43.94622     0
## 320  -69.72242 43.94622     0
## 321  -69.72242 43.94622     0
## 322  -69.72242 43.94622     0
## 323  -74.07045 40.05046    -1
## 324  -74.07045 40.05046    -1
## 325  -74.07045 40.05046    -1
## 326  -74.07045 40.05046    -1
## 327  -74.07045 40.05046    -1
## 328  -74.07045 40.05046    -1
## 329  -74.07045 40.05046    -1
## 330  -74.07045 40.05046    -1
## 331  -74.07045 40.05046    -1
## 332  -74.07045 40.05046    -1
## 333  -74.07045 40.05046    -1
## 334  -74.07045 40.05046    -1
## 335  -74.07045 40.05046    -1
## 336  -74.07045 40.05046    -1
## 337  -74.07045 40.05046    -1
## 338  -74.07045 40.05046    -1
## 339  -74.07045 40.05046    -1
## 340  -74.07045 40.05046    -1
## 341  -74.07045 40.05046    -1
## 342  -74.07045 40.05046    -1
## 343  -74.07045 40.05046    -1
## 344  -74.07045 40.05046    -1
## 345  -74.07045 40.05046    -1
## 346  -74.07045 40.05046    -1
## 347  -74.07045 40.05046    -1
## 348  -74.07045 40.05046    -1
## 349  -74.07045 40.05046    -1
## 350  -74.07045 40.05046    -1
## 351  -74.07045 40.05046    -1
## 352  -74.07045 40.05046    -1
## 353  -74.07045 40.05046    -1
## 354  -74.07045 40.05046    -1
## 355  -74.07045 40.05046    -1
## 356  -74.07045 40.05046    -1
## 357  -74.07045 40.05046    -1
## 358  -74.07045 40.05046    -1
## 359  -74.07045 40.05046    -1
## 360  -74.07045 40.05046    -1
## 361  -74.07045 40.05046    -1
## 362  -74.07045 40.05046    -1
## 363  -74.07045 40.05046    -1
## 364  -74.07045 40.05046    -1
## 365  -74.07045 40.05046    -1
## 366  -74.07045 40.05046    -1
## 367  -74.07045 40.05046    -1
## 368  -74.07045 40.05046    -1
## 369  -74.07045 40.05046    -1
## 370  -74.07045 40.05046    -1
## 371  -74.07045 40.05046    -1
## 372  -74.07045 40.05046    -1
## 373  -74.07045 40.05046    -1
## 374  -74.07045 40.05046    -1
## 375  -74.07045 40.05046    -1
## 376  -74.07045 40.05046    -1
## 377  -74.07045 40.05046    -1
## 378  -74.07045 40.05046    -1
## 379  -74.07045 40.05046    -1
## 380  -74.07045 40.05046    -1
## 381  -74.07045 40.05046    -1
## 382  -74.07045 40.05046    -1
## 383  -74.07045 40.05046    -1
## 384  -74.07045 40.05046    -1
## 385  -74.07045 40.05046    -1
## 386  -74.07045 40.05046    -1
## 387  -74.07045 40.05046    -1
## 388  -74.07045 40.05046    -1
## 389  -74.07045 40.05046    -1
## 390  -74.07045 40.05046    -1
## 391  -74.07045 40.05046    -1
## 392  -74.07045 40.05046    -1
## 393  -74.07045 40.05046    -1
## 394  -74.07045 40.05046    -1
## 395  -74.07045 40.05046    -1
## 396  -74.07045 40.05046    -1
## 397  -74.07045 40.05046    -1
## 398  -74.07045 40.05046    -1
## 399  -74.07045 40.05046    -1
## 400  -74.07045 40.05046    -1
## 401  -74.07045 40.05046    -1
## 402  -74.07045 40.05046    -1
## 403  -74.07045 40.05046    -1
## 404  -74.07045 40.05046    -1
## 405  -74.07045 40.05046    -1
## 406  -74.07045 40.05046    -1
## 407  -74.07045 40.05046    -1
## 408  -74.07045 40.05046    -1
## 409  -74.07045 40.05046    -1
## 410  -74.07045 40.05046    -1
## 411  -74.07045 40.05046    -1
## 412  -74.07045 40.05046    -1
## 413  -74.07045 40.05046    -1
## 414  -74.07045 40.05046    -1
## 415  -74.07045 40.05046    -1
## 416  -74.07045 40.05046    -1
## 417  -74.07045 40.05046    -1
## 418  -74.07045 40.05046    -1
## 419  -74.07045 40.05046    -1
## 420  -74.07045 40.05046    -1
## 421  -74.07045 40.05046    -1
## 422  -74.07045 40.05046    -1
## 423  -74.07045 40.05046    -1
## 424  -74.07045 40.05046    -1
## 425  -74.07045 40.05046    -1
## 426  -74.07045 40.05046    -1
## 427  -74.07045 40.05046    -1
## 428  -74.07045 40.05046    -1
## 429  -74.07045 40.05046    -1
## 430  -74.07045 40.05046    -1
## 431  -74.07045 40.05046    -1
## 432  -74.07045 40.05046    -1
## 433  -74.07045 40.05046    -1
## 434  -74.07045 40.05046    -1
## 435  -74.07045 40.05046    -1
## 436  -74.07045 40.05046    -1
## 437  -74.07045 40.05046    -1
## 438  -74.07045 40.05046    -1
## 439  -74.07045 40.05046    -1
## 440  -74.07045 40.05046    -1
## 441  -74.07045 40.05046    -1
## 442  -74.07045 40.05046    -1
## 443  -74.07045 40.05046    -1
## 444  -74.07045 40.05046    -1
## 445  -74.07045 40.05046    -1
## 446  -74.07045 40.05046    -1
## 447  -74.07045 40.05046    -1
## 448  -74.07045 40.05046    -1
## 449  -74.07045 40.05046    -1
## 450  -74.07045 40.05046    -1
## 451  -74.07045 40.05046    -1
## 452  -74.07045 40.05046    -1
## 453  -74.07045 40.05046    -1
## 454  -74.07045 40.05046    -1
## 455  -74.07045 40.05046    -1
## 456  -74.07045 40.05046    -1
## 457  -74.07045 40.05046    -1
## 458  -74.07045 40.05046    -1
## 459  -74.07045 40.05046    -1
## 460  -74.07045 40.05046    -1
## 461  -74.07045 40.05046    -1
## 462  -74.07045 40.05046    -1
## 463  -74.07045 40.05046    -1
## 464  -74.07045 40.05046    -1
## 465  -74.07045 40.05046    -1
## 466  -74.07045 40.05046    -1
## 467  -74.07045 40.05046    -1
## 468  -74.07045 40.05046    -1
## 469  -74.07045 40.05046    -1
## 470  -74.07045 40.05046    -1
## 471  -74.07045 40.05046    -1
## 472  -74.07045 40.05046    -1
## 473  -74.07045 40.05046    -1
## 474  -74.07045 40.05046    -1
## 475  -74.07045 40.05046    -1
## 476  -74.07045 40.05046    -1
## 477  -74.07045 40.05046    -1
## 478  -74.07045 40.05046    -1
## 479  -74.07045 40.05046    -1
## 480  -74.07045 40.05046    -1
## 481  -74.07045 40.05046    -1
## 482  -74.07045 40.05046    -1
## 483  -74.07045 40.05046    -1
## 484  -74.07045 40.05046    -1
## 485  -74.07045 40.05046    -1
## 486  -74.07045 40.05046    -1
## 487  -74.07045 40.05046    -1
## 488  -74.07045 40.05046    -1
## 489  -74.07045 40.05046    -1
## 490  -74.07045 40.05046    -1
## 491  -74.07045 40.05046    -1
## 492  -74.07045 40.05046    -1
## 493  -74.07045 40.05046    -1
## 494  -74.07045 40.05046    -1
## 495  -74.07045 40.05046    -1
## 496  -74.07045 40.05046    -1
## 497  -74.07045 40.05046    -1
## 498  -74.07045 40.05046    -1
## 499  -74.07045 40.05046    -1
## 500  -74.07045 40.05046    -1
## 501  -74.07045 40.05046    -1
## 502  -74.07045 40.05046    -1
## 503  -74.07045 40.05046    -1
## 504  -74.07045 40.05046    -1
## 505  -74.07045 40.05046    -1
## 506  -74.07045 40.05046    -1
## 507  -74.07045 40.05046    -1
## 508  -74.07045 40.05046    -1
## 509  -74.07045 40.05046    -1
## 510  -74.07045 40.05046    -1
## 511  -74.07045 40.05046    -1
## 512  -74.07045 40.05046    -1
## 513  -74.07045 40.05046    -1
## 514  -74.07045 40.05046    -1
## 515  -74.07045 40.05046    -1
## 516  -74.07045 40.05046    -1
## 517  -74.07045 40.05046    -1
## 518  -74.07045 40.05046    -1
## 519  -74.07045 40.05046    -1
## 520  -74.07045 40.05046    -1
## 521  -74.07045 40.05046    -1
## 522  -74.07045 40.05046    -1
## 523  -74.07045 40.05046    -1
## 524  -74.07045 40.05046    -1
## 525  -74.07045 40.05046    -1
## 526  -74.07045 40.05046    -1
## 527  -74.07045 40.05046    -1
## 528  -74.07045 40.05046    -1
## 529  -74.07045 40.05046    -1
## 530  -74.07045 40.05046    -1
## 531  -74.07045 40.05046    -1
## 532  -74.07045 40.05046    -1
## 533  -74.07045 40.05046    -1
## 534  -74.07045 40.05046    -1
## 535  -74.07045 40.05046    -1
## 536  -74.07045 40.05046    -1
## 537  -74.07045 40.05046    -1
## 538  -74.07045 40.05046    -1
## 539  -74.07045 40.05046    -1
## 540  -74.07045 40.05046    -1
## 541  -74.07045 40.05046    -1
## 542  -74.07045 40.05046    -1
## 543  -74.07045 40.05046    -1
## 544  -74.07045 40.05046    -1
## 545  -74.07045 40.05046    -1
## 546  -74.07045 40.05046    -1
## 547  -74.07045 40.05046    -1
## 548  -74.07045 40.05046    -1
## 549  -74.07045 40.05046    -1
## 550  -74.07045 40.05046    -1
## 551  -74.07045 40.05046    -1
## 552  -74.07045 40.05046    -1
## 553  -74.07045 40.05046    -1
## 554  -74.07045 40.05046    -1
## 555  -74.07045 40.05046    -1
## 556  -74.07045 40.05046    -1
## 557  -74.07045 40.05046    -1
## 558  -74.07045 40.05046    -1
## 559  -74.07045 40.05046    -1
## 560  -74.07045 40.05046    -1
## 561  -74.07045 40.05046    -1
## 562  -74.07045 40.05046    -1
## 563  -74.07045 40.05046    -1
## 564  -74.07045 40.05046    -1
## 565  -74.07045 40.05046    -1
## 566  -74.07045 40.05046    -1
## 567  -74.07045 40.05046    -1
## 568  -74.07045 40.05046    -1
## 569  -74.07045 40.05046    -1
## 570  -74.07045 40.05046    -1
## 571  -74.07045 40.05046    -1
## 572  -74.07045 40.05046    -1
## 573  -74.07045 40.05046    -1
## 574  -74.07045 40.05046    -1
## 575  -74.07045 40.05046    -1
## 576  -74.07045 40.05046    -1
## 577  -74.07045 40.05046    -1
## 578  -74.07045 40.05046    -1
## 579  -74.07045 40.05046    -1
## 580  -74.07045 40.05046    -1
## 581  -74.07045 40.05046    -1
## 582  -74.07045 40.05046    -1
## 583  -74.07045 40.05046    -1
## 584  -74.07045 40.05046    -1
## 585  -74.07045 40.05046    -1
## 586  -74.07045 40.05046    -1
## 587  -74.07045 40.05046    -1
## 588  -74.07045 40.05046    -1
## 589  -74.07045 40.05046    -1
## 590  -74.07045 40.05046    -1
## 591  -74.07045 40.05046    -1
## 592  -74.07045 40.05046    -1
## 593  -74.07045 40.05046    -1
## 594  -74.07045 40.05046    -1
## 595  -74.07045 40.05046    -1
## 596  -74.07045 40.05046    -1
## 597  -74.07045 40.05046    -1
## 598  -74.07045 40.05046    -1
## 599  -74.07045 40.05046    -1
## 600  -74.07045 40.05046    -1
## 601  -74.07045 40.05046    -1
## 602  -74.07045 40.05046    -1
## 603  -74.07045 40.05046    -1
## 604  -74.07045 40.05046    -1
## 605  -74.07045 40.05046    -1
## 606  -74.07045 40.05046    -1
## 607  -74.07045 40.05046    -1
## 608  -74.07045 40.05046    -1
## 609  -74.07045 40.05046    -1
## 610  -74.07045 40.05046    -1
## 611  -74.07045 40.05046    -1
## 612  -74.07045 40.05046    -1
## 613  -74.07045 40.05046    -1
## 614  -74.07045 40.05046    -1
## 615  -74.07045 40.05046    -1
## 616  -74.07045 40.05046    -1
## 617  -74.07045 40.05046    -1
## 618  -74.07045 40.05046    -1
## 619  -74.07045 40.05046    -1
## 620  -74.07045 40.05046    -1
## 621  -74.07045 40.05046    -1
## 622  -74.07045 40.05046    -1
## 623  -74.07045 40.05046    -1
## 624  -74.07045 40.05046    -1
## 625  -74.07045 40.05046    -1
## 626  -74.07045 40.05046    -1
## 627  -74.07045 40.05046    -1
## 628  -74.07045 40.05046    -1
## 629  -74.07045 40.05046    -1
## 630  -74.07045 40.05046    -1
## 631  -74.07045 40.05046    -1
## 632  -74.07045 40.05046    -1
## 633  -74.07045 40.05046    -1
## 634  -74.07045 40.05046    -1
## 635  -74.07045 40.05046    -1
## 636  -74.07045 40.05046    -1
## 637  -70.90760 41.62486    -6
## 638  -70.90760 41.62486    -6
## 639  -70.90760 41.62486    -6
## 640  -70.90760 41.62486    -6
## 641  -70.90760 41.62486    -6
## 642  -70.90760 41.62486    -6
## 643  -70.90760 41.62486    -6
## 644  -70.90760 41.62486    -6
## 645  -70.90760 41.62486    -6
## 646  -70.90760 41.62486    -6
## 647  -70.90760 41.62486    -6
## 648  -70.90760 41.62486    -6
## 649  -70.90760 41.62486    -6
## 650  -70.90760 41.62486    -6
## 651  -70.90760 41.62486    -6
## 652  -70.90760 41.62486    -6
## 653  -70.90760 41.62486    -6
## 654  -70.90760 41.62486    -6
## 655  -70.90760 41.62486    -6
## 656  -70.90760 41.62486    -6
## 657  -70.90760 41.62486    -6
## 658  -70.90760 41.62486    -6
## 659  -70.90760 41.62486    -6
## 660  -70.90760 41.62486    -6
## 661  -70.90760 41.62486    -6
## 662  -70.90760 41.62486    -6
## 663  -70.90760 41.62486    -6
## 664  -70.90760 41.62486    -6
## 665  -70.90760 41.62486    -6
## 666  -70.90760 41.62486    -6
## 667  -75.65578 35.91287     0
## 668  -75.65578 35.91287     0
## 669  -75.65578 35.91287     0
## 670  -75.65578 35.91287     0
## 671  -75.65578 35.91287     0
## 672  -75.65578 35.91287     0
## 673  -75.65578 35.91287     0
## 674  -75.65578 35.91287     0
## 675  -75.65578 35.91287     0
## 676  -75.65578 35.91287     0
## 677  -75.65578 35.91287     0
## 678  -75.65578 35.91287     0
## 679  -75.65578 35.91287     0
## 680  -75.65578 35.91287     0
## 681  -75.65578 35.91287     0
## 682  -75.65578 35.91287     0
## 683  -75.65578 35.91287     0
## 684  -75.65578 35.91287     0
## 685  -75.65578 35.91287     0
## 686  -75.65578 35.91287     0
## 687  -75.65578 35.91287     0
## 688  -75.65578 35.91287     0
## 689  -74.18475 39.80866     1
## 690  -74.18475 39.80866     1
## 691  -74.18475 39.80866     1
## 692  -74.18475 39.80866     1
## 693  -74.18475 39.80866     1
## 694  -74.18475 39.80866     1
## 695  -74.18475 39.80866     1
## 696  -74.18475 39.80866     1
## 697  -74.18475 39.80866     1
## 698  -74.18475 39.80866     1
## 699  -74.18475 39.80866     1
## 700  -74.18475 39.80866     1
## 701  -74.18475 39.80866     1
## 702  -74.18475 39.80866     1
## 703  -74.18475 39.80866     1
## 704  -74.18475 39.80866     1
## 705  -74.18475 39.80866     1
## 706  -74.18475 39.80866     1
## 707  -74.18475 39.80866     1
## 708  -74.18475 39.80866     1
## 709  -74.18475 39.80866     1
## 710  -74.18475 39.80866     1
## 711  -74.18475 39.80866     1
## 712  -74.18475 39.80866     1
## 713  -74.18475 39.80866     1
## 714  -74.18475 39.80866     1
## 715  -74.18475 39.80866     1
## 716  -74.18475 39.80866     1
## 717  -74.18475 39.80866     1
## 718  -74.18475 39.80866     1
## 719  -74.18475 39.80866     1
## 720  -74.18475 39.80866     1
## 721  -74.18475 39.80866     1
## 722  -74.18475 39.80866     1
## 723  -74.18475 39.80866     1
## 724  -74.18475 39.80866     1
## 725  -74.18475 39.80866     1
## 726  -74.18475 39.80866     1
## 727  -74.18475 39.80866     1
## 728  -74.18475 39.80866     1
## 729  -74.18475 39.80866     1
## 730  -74.18475 39.80866     1
## 731  -74.18475 39.80866     1
## 732  -70.90760 41.62486    -6
## 733  -70.90760 41.62486    -6
## 734  -70.90760 41.62486    -6
## 735  -70.90760 41.62486    -6
## 736  -70.90760 41.62486    -6
## 737  -70.90760 41.62486    -6
## 738  -70.90760 41.62486    -6
## 739  -70.90760 41.62486    -6
## 740  -70.90760 41.62486    -6
## 741  -70.90760 41.62486    -6
## 742  -70.90760 41.62486    -6
## 743  -70.90760 41.62486    -6
## 744  -70.90760 41.62486    -6
## 745  -70.90760 41.62486    -6
## 746  -70.90760 41.62486    -6
## 747  -70.90760 41.62486    -6
## 748  -70.90760 41.62486    -6
## 749  -70.90760 41.62486    -6
## 750  -70.90760 41.62486    -6
## 751  -70.90760 41.62486    -6
## 752  -70.90760 41.62486    -6
## 753  -70.90760 41.62486    -6
## 754  -70.90760 41.62486    -6
## 755  -70.90760 41.62486    -6
## 756  -70.90760 41.62486    -6
## 757  -70.90760 41.62486    -6
## 758  -70.90760 41.62486    -6
## 759  -70.90760 41.62486    -6
## 760  -70.90760 41.62486    -6
## 761  -70.90760 41.62486    -6
## 762  -74.32448 39.50878    -3
## 763  -74.32448 39.50878    -3
## 764  -74.32448 39.50878    -3
## 765  -74.32448 39.50878    -3
## 766  -74.32448 39.50878    -3
## 767  -74.32448 39.50878    -3
## 768  -74.32448 39.50878    -3
## 769  -74.32448 39.50878    -3
## 770  -74.32448 39.50878    -3
## 771  -74.32448 39.50878    -3
## 772  -74.32448 39.50878    -3
## 773  -74.32448 39.50878    -3
## 774  -74.32448 39.50878    -3
## 775  -74.32448 39.50878    -3
## 776  -74.32448 39.50878    -3
## 777  -74.32448 39.50878    -3
## 778  -74.32448 39.50878    -3
## 779  -74.32448 39.50878    -3
## 780  -74.32448 39.50878    -3
## 781  -74.32448 39.50878    -3
## 782  -74.32448 39.50878    -3
## 783  -74.32448 39.50878    -3
## 784  -74.32448 39.50878    -3
## 785  -74.32448 39.50878    -3
## 786  -74.32448 39.50878    -3
## 787  -74.32448 39.50878    -3
## 788  -74.32448 39.50878    -3
## 789  -74.32448 39.50878    -3
## 790  -74.32448 39.50878    -3
## 791  -74.32448 39.50878    -3
## 792  -74.32448 39.50878    -3
## 793  -74.32448 39.50878    -3
## 794  -74.32448 39.50878    -3
## 795  -74.32448 39.50878    -3
## 796  -74.32448 39.50878    -3
## 797  -74.32448 39.50878    -3
## 798  -74.32448 39.50878    -3
## 799  -74.32448 39.50878    -3
## 800  -74.32448 39.50878    -3
## 801  -74.32448 39.50878    -3
## 802  -74.32448 39.50878    -3
## 803  -74.32448 39.50878    -3
## 804  -74.32448 39.50878    -3
## 805  -74.32448 39.50878    -3
## 806  -74.32448 39.50878    -3
## 807  -74.32448 39.50878    -3
## 808  -74.32448 39.50878    -3
## 809  -74.32448 39.50878    -3
## 810  -74.32448 39.50878    -3
## 811  -74.32448 39.50878    -3
## 812  -74.32448 39.50878    -3
## 813  -74.32448 39.50878    -3
## 814  -74.32448 39.50878    -3
## 815  -74.32448 39.50878    -3
## 816  -74.32448 39.50878    -3
## 817  -74.32448 39.50878    -3
## 818  -74.32448 39.50878    -3
## 819  -74.32448 39.50878    -3
## 820  -74.32448 39.50878    -3
## 821  -74.32448 39.50878    -3
## 822  -74.32448 39.50878    -3
## 823  -74.32448 39.50878    -3
## 824  -74.32448 39.50878    -3
## 825  -74.32448 39.50878    -3
## 826  -74.32448 39.50878    -3
## 827  -74.32448 39.50878    -3
## 828  -74.32448 39.50878    -3
## 829  -74.32448 39.50878    -3
## 830  -74.32448 39.50878    -3
## 831  -74.32448 39.50878    -3
## 832  -74.32448 39.50878    -3
## 833  -74.32448 39.50878    -3
## 834  -74.32448 39.50878    -3
## 835  -74.32448 39.50878    -3
## 836  -74.32448 39.50878    -3
## 837  -74.32448 39.50878    -3
## 838  -74.32448 39.50878    -3
## 839  -74.32448 39.50878    -3
## 840  -74.32448 39.50878    -3
## 841  -74.32448 39.50878    -3
## 842  -74.32448 39.50878    -3
## 843  -74.32448 39.50878    -3
## 844  -74.32448 39.50878    -3
## 845  -74.32448 39.50878    -3
## 846  -74.32448 39.50878    -3
## 847  -74.32448 39.50878    -3
## 848  -74.32448 39.50878    -3
## 849  -74.32448 39.50878    -3
## 850  -74.32448 39.50878    -3
## 851  -74.32448 39.50878    -3
## 852  -74.32448 39.50878    -3
## 853  -74.32448 39.50878    -3
## 854  -74.32448 39.50878    -3
## 855  -74.32448 39.50878    -3
## 856  -74.32448 39.50878    -3
## 857  -74.32448 39.50878    -3
## 858  -74.32448 39.50878    -3
## 859  -74.32448 39.50878    -3
## 860  -74.32448 39.50878    -3
## 861  -74.32448 39.50878    -3
## 862  -74.32448 39.50878    -3
## 863  -74.32448 39.50878    -3
## 864  -74.32448 39.50878    -3
## 865  -74.32448 39.50878    -3
## 866  -74.32448 39.50878    -3
## 867  -74.32448 39.50878    -3
## 868  -74.32448 39.50878    -3
## 869  -74.32448 39.50878    -3
## 870  -74.32448 39.50878    -3
## 871  -74.32448 39.50878    -3
## 872  -74.32448 39.50878    -3
## 873  -74.32448 39.50878    -3
## 874  -74.32448 39.50878    -3
## 875  -74.32448 39.50878    -3
## 876  -74.32448 39.50878    -3
## 877  -74.32448 39.50878    -3
## 878  -74.32448 39.50878    -3
## 879  -74.32448 39.50878    -3
## 880  -74.32448 39.50878    -3
## 881  -74.32448 39.50878    -3
## 882  -74.32448 39.50878    -3
## 883  -74.32448 39.50878    -3
## 884  -74.32448 39.50878    -3
## 885  -74.32448 39.50878    -3
## 886  -74.32448 39.50878    -3
## 887  -74.32448 39.50878    -3
## 888  -74.32448 39.50878    -3
## 889  -74.32448 39.50878    -3
## 890  -74.32448 39.50878    -3
## 891  -74.32448 39.50878    -3
## 892  -74.32448 39.50878    -3
## 893  -74.32448 39.50878    -3
## 894  -74.32448 39.50878    -3
## 895  -74.32448 39.50878    -3
## 896  -74.32448 39.50878    -3
## 897  -74.32448 39.50878    -3
## 898  -74.32448 39.50878    -3
## 899  -74.32448 39.50878    -3
## 900  -74.32448 39.50878    -3
## 901  -74.32448 39.50878    -3
## 902  -74.32448 39.50878    -3
## 903  -74.32448 39.50878    -3
## 904  -74.32448 39.50878    -3
## 905  -74.32448 39.50878    -3
## 906  -74.32448 39.50878    -3
## 907  -74.32448 39.50878    -3
## 908  -74.32448 39.50878    -3
## 909  -74.32448 39.50878    -3
## 910  -74.32448 39.50878    -3
## 911  -74.32448 39.50878    -3
## 912  -74.32448 39.50878    -3
## 913  -74.32448 39.50878    -3
## 914  -74.32448 39.50878    -3
## 915  -74.32448 39.50878    -3
## 916  -74.32448 39.50878    -3
## 917  -74.32448 39.50878    -3
## 918  -74.32448 39.50878    -3
## 919  -74.32448 39.50878    -3
## 920  -74.32448 39.50878    -3
## 921  -74.32448 39.50878    -3
## 922  -74.32448 39.50878    -3
## 923  -74.32448 39.50878    -3
## 924  -74.32448 39.50878    -3
## 925  -74.32448 39.50878    -3
## 926  -74.32448 39.50878    -3
## 927  -74.32448 39.50878    -3
## 928  -74.32448 39.50878    -3
## 929  -74.32448 39.50878    -3
## 930  -74.32448 39.50878    -3
## 931  -74.32448 39.50878    -3
## 932  -74.32448 39.50878    -3
## 933  -74.32448 39.50878    -3
## 934  -74.32448 39.50878    -3
## 935  -74.32448 39.50878    -3
## 936  -74.32448 39.50878    -3
## 937  -74.32448 39.50878    -3
## 938  -74.32448 39.50878    -3
## 939  -74.32448 39.50878    -3
## 940  -74.32448 39.50878    -3
## 941  -74.32448 39.50878    -3
## 942  -74.32448 39.50878    -3
## 943  -74.32448 39.50878    -3
## 944  -74.32448 39.50878    -3
## 945  -74.32448 39.50878    -3
## 946  -74.32448 39.50878    -3
## 947  -74.32448 39.50878    -3
## 948  -74.32448 39.50878    -3
## 949  -74.32448 39.50878    -3
## 950  -74.32448 39.50878    -3
## 951  -74.32448 39.50878    -3
## 952  -74.32448 39.50878    -3
## 953  -74.32448 39.50878    -3
## 954  -74.32448 39.50878    -3
## 955  -74.32448 39.50878    -3
## 956  -74.32448 39.50878    -3
## 957  -74.32448 39.50878    -3
## 958  -74.32448 39.50878    -3
## 959  -74.32448 39.50878    -3
## 960  -74.32448 39.50878    -3
## 961  -74.32448 39.50878    -3
## 962  -74.32448 39.50878    -3
## 963  -74.32448 39.50878    -3
## 964  -74.32448 39.50878    -3
## 965  -74.32448 39.50878    -3
## 966  -74.32448 39.50878    -3
## 967  -74.32448 39.50878    -3
## 968  -74.32448 39.50878    -3
## 969  -74.32448 39.50878    -3
## 970  -74.32448 39.50878    -3
## 971  -74.32448 39.50878    -3
## 972  -74.32448 39.50878    -3
## 973  -74.32448 39.50878    -3
## 974  -74.32448 39.50878    -3
## 975  -74.32448 39.50878    -3
## 976  -74.32448 39.50878    -3
## 977  -74.32448 39.50878    -3
## 978  -74.32448 39.50878    -3
## 979  -74.32448 39.50878    -3
## 980  -74.32448 39.50878    -3
## 981  -74.32448 39.50878    -3
## 982  -74.32448 39.50878    -3
## 983  -74.32448 39.50878    -3
## 984  -74.32448 39.50878    -3
## 985  -74.32448 39.50878    -3
## 986  -74.32448 39.50878    -3
## 987  -74.32448 39.50878    -3
## 988  -74.32448 39.50878    -3
## 989  -74.32448 39.50878    -3
## 990  -74.32448 39.50878    -3
## 991  -74.32448 39.50878    -3
## 992  -74.32448 39.50878    -3
## 993  -74.32448 39.50878    -3
## 994  -74.32448 39.50878    -3
## 995  -74.32448 39.50878    -3
## 996  -74.32448 39.50878    -3
## 997  -74.32448 39.50878    -3
## 998  -74.32448 39.50878    -3
## 999  -74.32448 39.50878    -3
## 1000 -74.32448 39.50878    -3
## 1001 -74.32448 39.50878    -3
## 1002 -74.32448 39.50878    -3
## 1003 -74.32448 39.50878    -3
## 1004 -74.32448 39.50878    -3
## 1005 -74.32448 39.50878    -3
## 1006 -74.32448 39.50878    -3
## 1007 -74.32448 39.50878    -3
## 1008 -74.32448 39.50878    -3
## 1009 -74.32448 39.50878    -3
## 1010 -74.32448 39.50878    -3
## 1011 -74.32448 39.50878    -3
## 1012 -74.32448 39.50878    -3
## 1013 -74.32448 39.50878    -3
## 1014 -74.32448 39.50878    -3
## 1015 -74.32448 39.50878    -3
## 1016 -74.32448 39.50878    -3
## 1017 -74.32448 39.50878    -3
## 1018 -74.32448 39.50878    -3
## 1019 -74.32448 39.50878    -3
## 1020 -74.32448 39.50878    -3
## 1021 -74.32448 39.50878    -3
## 1022 -74.32448 39.50878    -3
## 1023 -74.32448 39.50878    -3
## 1024 -74.32448 39.50878    -3
## 1025 -74.32448 39.50878    -3
## 1026 -74.32448 39.50878    -3
## 1027 -74.32448 39.50878    -3
## 1028 -74.32448 39.50878    -3
## 1029 -74.32448 39.50878    -3
## 1030 -74.32448 39.50878    -3
## 1031 -74.32448 39.50878    -3
## 1032 -74.32448 39.50878    -3
## 1033 -74.32448 39.50878    -3
## 1034 -74.32448 39.50878    -3
## 1035 -74.32448 39.50878    -3
## 1036 -74.32448 39.50878    -3
## 1037 -74.32448 39.50878    -3
## 1038 -74.32448 39.50878    -3
## 1039 -74.32448 39.50878    -3
## 1040 -74.32448 39.50878    -3
## 1041 -74.32448 39.50878    -3
## 1042 -74.32448 39.50878    -3
## 1043 -74.32448 39.50878    -3
## 1044 -74.32448 39.50878    -3
## 1045 -74.32448 39.50878    -3
## 1046 -74.32448 39.50878    -3
## 1047 -74.32448 39.50878    -3
## 1048 -74.32448 39.50878    -3
## 1049 -74.32448 39.50878    -3
## 1050 -74.32448 39.50878    -3
## 1051 -74.32448 39.50878    -3
## 1052 -74.32448 39.50878    -3
## 1053 -74.32448 39.50878    -3
## 1054 -74.32448 39.50878    -3
## 1055 -74.32448 39.50878    -3
## 1056 -74.32448 39.50878    -3
## 1057 -74.32448 39.50878    -3
## 1058 -74.32448 39.50878    -3
## 1059 -74.32448 39.50878    -3
## 1060 -74.32448 39.50878    -3
## 1061 -74.32448 39.50878    -3
## 1062 -74.32448 39.50878    -3
## 1063 -74.32448 39.50878    -3
## 1064 -74.32448 39.50878    -3
## 1065 -74.32448 39.50878    -3
## 1066 -74.32448 39.50878    -3
## 1067 -74.32448 39.50878    -3
## 1068 -74.32448 39.50878    -3
## 1069 -74.32448 39.50878    -3
## 1070 -74.32448 39.50878    -3
## 1071 -74.32448 39.50878    -3
## 1072 -74.32448 39.50878    -3
## 1073 -74.32448 39.50878    -3
## 1074 -74.32448 39.50878    -3
## 1075 -74.32448 39.50878    -3
## 1076 -74.32448 39.50878    -3
## 1077 -74.32448 39.50878    -3
## 1078 -74.32448 39.50878    -3
## 1079 -74.32448 39.50878    -3
## 1080 -74.32448 39.50878    -3
## 1081 -74.32448 39.50878    -3
## 1082 -74.32448 39.50878    -3
## 1083 -74.32448 39.50878    -3
## 1084 -74.32448 39.50878    -3
## 1085 -74.32448 39.50878    -3
## 1086 -74.32448 39.50878    -3
## 1087 -74.32448 39.50878    -3
## 1088 -74.32448 39.50878    -3
## 1089 -74.32448 39.50878    -3
## 1090 -74.32448 39.50878    -3
## 1091 -74.32448 39.50878    -3
## 1092 -74.32448 39.50878    -3
## 1093 -74.32448 39.50878    -3
## 1094 -74.32448 39.50878    -3
## 1095 -74.32448 39.50878    -3
## 1096 -74.32448 39.50878    -3
## 1097 -74.32448 39.50878    -3
## 1098 -74.32448 39.50878    -3
## 1099 -74.32448 39.50878    -3
## 1100 -74.32448 39.50878    -3
## 1101 -74.32448 39.50878    -3
## 1102 -74.32448 39.50878    -3
## 1103 -74.32448 39.50878    -3
## 1104 -74.32448 39.50878    -3
## 1105 -74.32448 39.50878    -3
## 1106 -74.32448 39.50878    -3
## 1107 -74.32448 39.50878    -3
## 1108 -74.32448 39.50878    -3
## 1109 -74.32448 39.50878    -3
## 1110 -74.32448 39.50878    -3
## 1111 -74.32448 39.50878    -3
## 1112 -74.32448 39.50878    -3
## 1113 -74.32448 39.50878    -3
## 1114 -74.32448 39.50878    -3
## 1115 -74.32448 39.50878    -3
## 1116 -74.32448 39.50878    -3
## 1117 -74.32448 39.50878    -3
## 1118 -74.32448 39.50878    -3
## 1119 -74.32448 39.50878    -3
## 1120 -74.32448 39.50878    -3
## 1121 -74.32448 39.50878    -3
## 1122 -74.32448 39.50878    -3
## 1123 -74.32448 39.50878    -3
## 1124 -74.32448 39.50878    -3
## 1125 -74.32448 39.50878    -3
## 1126 -74.32448 39.50878    -3
## 1127 -74.32448 39.50878    -3
## 1128 -74.32448 39.50878    -3
## 1129 -74.32448 39.50878    -3
## 1130 -74.32448 39.50878    -3
## 1131 -74.32448 39.50878    -3
## 1132 -74.32448 39.50878    -3
## 1133 -74.32448 39.50878    -3
## 1134 -74.32448 39.50878    -3
## 1135 -74.32448 39.50878    -3
## 1136 -74.32448 39.50878    -3
## 1137 -74.32448 39.50878    -3
## 1138 -74.32448 39.50878    -3
## 1139 -74.32448 39.50878    -3
## 1140 -74.32448 39.50878    -3
## 1141 -74.32448 39.50878    -3
## 1142 -74.32448 39.50878    -3
## 1143 -74.32448 39.50878    -3
## 1144 -74.32448 39.50878    -3
## 1145 -74.32448 39.50878    -3
## 1146 -74.32448 39.50878    -3
## 1147 -74.32448 39.50878    -3
## 1148 -74.32448 39.50878    -3
## 1149 -74.32448 39.50878    -3
## 1150 -74.32448 39.50878    -3
## 1151 -74.32448 39.50878    -3
## 1152 -74.32448 39.50878    -3
## 1153 -74.32448 39.50878    -3
## 1154 -79.91624 32.75873    -5
## 1155 -79.91624 32.75873    -5
## 1156 -79.91624 32.75873    -5
## 1157 -79.91624 32.75873    -5
## 1158 -79.91624 32.75873    -5
## 1159 -79.91624 32.75873    -5
## 1160 -79.91624 32.75873    -5
## 1161 -79.91624 32.75873    -5
## 1162 -79.91624 32.75873    -5
## 1163 -79.91624 32.75873    -5
## 1164 -79.91624 32.75873    -5
## 1165 -79.91624 32.75873    -5
## 1166 -79.91624 32.75873    -5
## 1167 -79.91624 32.75873    -5
## 1168 -79.91624 32.75873    -5
## 1169 -79.91624 32.75873    -5
## 1170 -79.91624 32.75873    -5
## 1171 -79.91624 32.75873    -5
## 1172 -79.91624 32.75873    -5
## 1173 -79.91624 32.75873    -5
## 1174 -79.91624 32.75873    -5
## 1175 -79.91624 32.75873    -5
## 1176 -79.91624 32.75873    -5
## 1177 -74.75746 39.06786    -3
## 1178 -74.75746 39.06786    -3
## 1179 -74.75746 39.06786    -3
## 1180 -74.75746 39.06786    -3
## 1181 -74.75746 39.06786    -3
## 1182 -74.75746 39.06786    -3
## 1183 -74.75746 39.06786    -3
## 1184 -74.75746 39.06786    -3
## 1185 -74.75746 39.06786    -3
## 1186 -74.75746 39.06786    -3
## 1187 -74.75746 39.06786    -3
## 1188 -74.75746 39.06786    -3
## 1189 -74.75746 39.06786    -3
## 1190 -74.75746 39.06786    -3
## 1191 -74.75746 39.06786    -3
## 1192 -74.75746 39.06786    -3
## 1193 -74.75746 39.06786    -3
## 1194 -74.75746 39.06786    -3
## 1195 -74.75746 39.06786    -3
## 1196 -74.75746 39.06786    -3
## 1197 -74.75746 39.06786    -3
## 1198 -74.75746 39.06786    -3
## 1199 -74.75746 39.06786    -3
## 1200 -74.75746 39.06786    -3
## 1201 -74.75746 39.06786    -3
## 1202 -74.75746 39.06786    -3
## 1203 -74.75746 39.06786    -3
## 1204 -74.75746 39.06786    -3
## 1205 -74.75746 39.06786    -3
## 1206 -74.75746 39.06786    -3
## 1207 -74.75746 39.06786    -3
## 1208 -74.75746 39.06786    -3
## 1209 -74.75746 39.06786    -3
## 1210 -74.75746 39.06786    -3
## 1211 -74.75746 39.06786    -3
## 1212 -74.75746 39.06786    -3
## 1213 -74.75746 39.06786    -3
## 1214 -74.75746 39.06786    -3
## 1215 -74.75746 39.06786    -3
## 1216 -74.75746 39.06786    -3
## 1217 -74.75746 39.06786    -3
## 1218 -74.75746 39.06786    -3
## 1219 -74.75746 39.06786    -3
## 1220 -74.75746 39.06786    -3
## 1221 -74.75746 39.06786    -3
## 1222 -74.75746 39.06786    -3
## 1223 -74.75746 39.06786    -3
## 1224 -74.75746 39.06786    -3
## 1225 -74.75746 39.06786    -3
## 1226 -74.75746 39.06786    -3
## 1227 -74.75746 39.06786    -3
## 1228 -74.75746 39.06786    -3
## 1229 -74.75746 39.06786    -3
## 1230 -74.75746 39.06786    -3
## 1231 -74.75746 39.06786    -3
## 1232 -74.75746 39.06786    -3
## 1233 -74.75746 39.06786    -3
## 1234 -74.75746 39.06786    -3
## 1235 -74.75746 39.06786    -3
## 1236 -70.97979 41.53583    -1
## 1237 -70.97979 41.53583    -1
## 1238 -70.97979 41.53583    -1
## 1239 -70.97979 41.53583    -1
## 1240 -70.97979 41.53583    -1
## 1241 -70.97979 41.53583    -1
## 1242 -70.97979 41.53583    -1
## 1243 -70.97979 41.53583    -1
## 1244 -70.97979 41.53583    -1
## 1245 -70.97979 41.53583    -1
## 1246 -70.97979 41.53583    -1
## 1247 -70.97979 41.53583    -1
## 1248 -70.97979 41.53583    -1
## 1249 -70.97979 41.53583    -1
## 1250 -70.97979 41.53583    -1
## 1251 -70.97979 41.53583    -1
## 1252 -70.97979 41.53583    -1
## 1253 -70.97979 41.53583    -1
## 1254 -70.97979 41.53583    -1
## 1255 -70.97979 41.53583    -1
## 1256 -70.97979 41.53583    -1
## 1257 -70.97979 41.53583    -1
## 1258 -70.97979 41.53583    -1
## 1259 -70.97979 41.53583    -1
## 1260 -70.97979 41.53583    -1
## 1261 -70.97979 41.53583    -1
## 1262 -70.97979 41.53583    -1
## 1263 -70.97979 41.53583    -1
## 1264 -70.97979 41.53583    -1
## 1265 -70.97979 41.53583    -1
## 1266 -71.52557 41.38235    -6
## 1267 -71.52557 41.38235    -6
## 1268 -71.52557 41.38235    -6
## 1269 -71.52557 41.38235    -6
## 1270 -71.52557 41.38235    -6
## 1271 -71.52557 41.38235    -6
## 1272 -71.52557 41.38235    -6
## 1273 -71.52557 41.38235    -6
## 1274 -71.52557 41.38235    -6
## 1275 -71.52557 41.38235    -6
## 1276 -71.52557 41.38235    -6
## 1277 -71.52557 41.38235    -6
## 1278 -71.52557 41.38235    -6
## 1279 -71.52557 41.38235    -6
## 1280 -71.52557 41.38235    -6
## 1281 -71.52557 41.38235    -6
## 1282 -71.52557 41.38235    -6
## 1283 -71.52557 41.38235    -6
## 1284 -71.52557 41.38235    -6
## 1285 -71.52557 41.38235    -6
## 1286 -71.52557 41.38235    -6
## 1287 -71.52557 41.38235    -6
## 1288 -71.52557 41.38235    -6
## 1289 -71.52557 41.38235    -6
## 1290 -71.52557 41.38235    -6
## 1291 -71.52557 41.38235    -6
## 1292 -71.52557 41.38235    -6
## 1293 -71.52557 41.38235    -6
## 1294 -71.52557 41.38235    -6
## 1295 -71.52557 41.38235    -6
## 1296 -71.52557 41.38235    -6
## 1297 -71.52557 41.38235    -6
## 1298 -71.52557 41.38235    -6
## 1299 -71.52557 41.38235    -6
## 1300 -71.52557 41.38235    -6
## 1301 -71.52557 41.38235    -6
## 1302 -71.52557 41.38235    -6
## 1303 -71.52557 41.38235    -6
## 1304 -71.52557 41.38235    -6
## 1305 -71.52557 41.38235    -6
## 1306 -71.52557 41.38235    -6
## 1307 -71.52557 41.38235    -6
## 1308 -71.52557 41.38235    -6
## 1309 -71.52557 41.38235    -6
## 1310 -71.52557 41.38235    -6
## 1311 -71.52557 41.38235    -6
## 1312 -71.52557 41.38235    -6
## 1313 -71.52557 41.38235    -6
## 1314 -70.90760 41.62486    -6
## 1315 -70.90760 41.62486    -6
## 1316 -70.90760 41.62486    -6
## 1317 -70.90760 41.62486    -6
## 1318 -70.90760 41.62486    -6
## 1319 -70.90760 41.62486    -6
## 1320 -70.90760 41.62486    -6
## 1321 -70.90760 41.62486    -6
## 1322 -70.90760 41.62486    -6
## 1323 -70.90760 41.62486    -6
## 1324 -70.90760 41.62486    -6
## 1325 -70.90760 41.62486    -6
## 1326 -70.90760 41.62486    -6
## 1327 -70.90760 41.62486    -6
## 1328 -70.90760 41.62486    -6
## 1329 -70.90760 41.62486    -6
## 1330 -70.90760 41.62486    -6
## 1331 -70.90760 41.62486    -6
## 1332 -70.90760 41.62486    -6
## 1333 -70.90760 41.62486    -6
## 1334 -70.90760 41.62486    -6
## 1335 -70.90760 41.62486    -6
## 1336 -70.90760 41.62486    -6
## 1337 -70.90760 41.62486    -6
## 1338 -70.90760 41.62486    -6
## 1339 -70.90760 41.62486    -6
## 1340 -70.90760 41.62486    -6
## 1341 -70.90760 41.62486    -6
## 1342 -70.90760 41.62486    -6
## 1343 -70.90760 41.62486    -6
## 1344 -77.92477 34.02835     0
## 1345 -77.92477 34.02835     0
## 1346 -77.92477 34.02835     0
## 1347 -77.92477 34.02835     0
## 1348 -77.92477 34.02835     0
## 1349 -77.92477 34.02835     0
## 1350 -77.92477 34.02835     0
## 1351 -77.92477 34.02835     0
## 1352 -77.92477 34.02835     0
## 1353 -77.92477 34.02835     0
## 1354 -77.92477 34.02835     0
## 1355 -77.92477 34.02835     0
## 1356 -77.92477 34.02835     0
## 1357 -77.92477 34.02835     0
## 1358 -77.92477 34.02835     0
## 1359 -77.92477 34.02835     0
## 1360 -77.92477 34.02835     0
## 1361 -77.92477 34.02835     0
## 1362 -77.92477 34.02835     0
## 1363 -77.92477 34.02835     0
## 1364 -77.92477 34.02835     0
## 1365 -77.92477 34.02835     0
## 1366 -77.92477 34.02835     0
## 1367 -77.92477 34.02835     0
## 1368 -77.92477 34.02835     0
## 1369 -77.92477 34.02835     0
## 1370 -77.92477 34.02835     0
## 1371 -77.92477 34.02835     0
## Warning in lc.dist(trans, ind_sites[, c(4, 3)], res = "dist"): One or more points are located outside of the [1;-15] depth range. You will get unrealistically huge distances unless you either increase the range of possible depths in trans.mat() or you move the problematic points in a spot where their depths fall within the [1;-15] depth range.
## You can use get.depth() to determine the depth of any point on a bathymetric map
#calculate dbmem and plot
dbmem_pos<-dbmem(distance_matrix, MEM.autocor = "positive")
## Warning in dudi.pco(matdist, scannf = FALSE, nf = 2): Non euclidean
## distance
plot(attributes(dbmem_pos)$values)
abline(h = 0)

The first 8 dbmems have vaules above 0 (i.e describe spatial autocorrelation)

temperatures

The environmental data is daily SSTs (NOAA OI SST V2 High Resolution Dataset) from

We extracted 10 years of daily optimum interpolated sea surface temperatures (OISSTs) from NOAA/OAR/ESRL PSD (Reynolds et al. 2007) for each of the sampling locations and calculated three derived environmental variables: mean temperature, mean annual minimum and mean annual maximum (table 5.1, fig. 5.1). The first principal component in a principal component analysis captured 90.1% of the variance among these variables, and all three variables were highly correlated as indicated by the biplot (supplementary fig. 5.1). While mean temperature and mean annual maximum demonstrate a linear relationship with latitude (R2 = 0.87 and 0.97; slope = -0.97 and - 128 1.1°C per ° latitude, respectively), minimum temperature is constrained by the freezing temperature of seawater and plateaued beyond 39 degrees of latitude. The value of the first temperature principal component and position along the shore for the sampling locations are also strongly correlated (R2 = 0.94, p < 10-9 ANOVA)

#env variable for pops
env<-read.table("./environmental/temp_summary.txt", sep = "\t", header = T)
env <- env %>%
  left_join(., pops, by = "pop")
## Warning: Column `pop` joining factor and character vector, coercing into
## character vector

genetic data

We use the covariance matrix calculated by PCANGSD to conduct a PCA that describes the genetic variation among samples. Then choose the number of PCs to retain using the Kaiser Guttman criteria

#pca
#to load pca: attach("./pcangsd/pca.r")
#snp_pca<-read.table("./pcangsd/pca.txt", header = T, sep= "\t")

#kaiser guttman
cutoff<-mean(pc1$sdev^2)
kg <- length((pc1$sdev^2)[(pc1$sdev^2)>cutoff])

#pcs to keep
snp_pcs <- pc1$x[,c(1:103)]

We keep the first kg pcs to fit our RDA on

RDA

Finally we get to running the RDA. We conduct a global RDA, an RDA for the temperature once conditioned on spatial data, and an RDA on the genetic DATA + spatial data once conditioned on the temperature. then we do model selection, significance testing and variance partitioning

#global rda
rda_null <- rda(snp_pcs ~ 1)
rda_ind<-rda(snp_pcs ~ as.matrix(env[,c(2)]) +as.matrix(env[,c(3)])  + as.matrix(env[,4]) +   (as.matrix(dbmem_pos[,1])) + (as.matrix(dbmem_pos[,2])) + (as.matrix(dbmem_pos[,3])) + (as.matrix(dbmem_pos[,4]))+ (as.matrix(dbmem_pos[,5])) + (as.matrix(dbmem_pos[,6])) + (as.matrix(dbmem_pos[,7])) + (as.matrix(dbmem_pos[,8])))

#variable selection
ordi <- ordistep(rda_null, scope = formula(rda_ind))
#ordi$anova

#variance partitioning
vp <- varpart(snp_pcs ,as.matrix(env[,c(2,3,4)]) , (as.matrix(dbmem_pos)[,c(1:8)]))

#controlling for spatial autocor
rda_ind_con_spatial<-rda(snp_pcs ~ as.matrix(env[,c(2)]) + Condition((as.matrix(dbmem_pos)[,c(1:8)])))

#controlling for temp
rda_ind_con_temp<-rda(snp_pcs ~ Condition(as.matrix(env[,c(2)])) + (as.matrix(dbmem_pos)[,c(1:8)]))

#anovas
aov_global <- anova(rda_ind)
aov_con_spatial <- anova(rda_ind_con_spatial) 
aov_con_temp <- anova(rda_ind_con_temp)

RDA Results Summary

This is a preliminary summary, more detailed results are to come.

We used the ordiR2step function to select our model. Using forward selection the full model was chosen (8 dbMEMs describing spatial autocorrelation and 3 environmental summary variables (mean temperature mean annual minimum, and mean annual maximum)).

Three RDAs were fit and tested for significance:
* global RDA - both the effect of spatial autocorrelation (the dbMEMs) and the effect of temerpature - p < 0.001
* RDA of temperature effect, conditioned on spatial autocorrelation - p < 0.001
* RDA of spatial effect conditioned on temperature - p < 0.001

Collectively these suggest that there is a significant relationship between the environmental temperature variable and the genetic covariance among samples, once the genetic data is conditioned on the effect of spatial autocorrelation AND vice versa. 76% of the total variation in the genetic dataset (or at least the first 103 of 1371 principal components of it) was constrained by the RDA.

Using the variance partitioning the relative effect of spatial autocorrelation and temperature could be estimated. Adjusted R2 for the effect of temperature, once spatial autocorrelation is accounted for, is 0.4%. In contrast 62% of the genetic dataset was explained by spatial autocorrelation alone. 14% of variation was explained jointly and could not be parsed. 23% was residual.

Examining the triplot suggests that the first dbmem separates northern from southern at the NJ.

RDA figures and results tables

#triplot
rda_sum<-summary(rda_ind)
scores<-as.data.frame(scores(rda_ind)$sites)
scores$pop<-pops$pop
arrows<-as.data.frame(rda_sum$biplot)
  #order the factors
scores$pop <- factor( as.character(scores$pop), levels=c("GA", "SC", "WNC", "NC", "KC", "SH", "RB", "OC", "MG", "CT", "SR", "BP", "HB",  "NBH", "P", "SLC", "ME", "MDIBL")  )
ggplot()+geom_point(aes(RDA1, RDA2, color = pop), data = scores)+geom_segment(data=arrows, aes(x=0, y=0, xend=RDA1, yend=RDA2), arrow=arrow(length=unit(0.2,"cm")), alpha=0.75, color="red")+geom_text(data=arrows, aes(x=RDA1+0.2, y=RDA2, label=c("Max Temp", "Min Temp", "Mean Temp", "dbMEM1", "dbMEM2","dbMEM3","dbMEM4","dbMEM5", "dbMEM6","dbMEM7", "dbMEM8")), size = 5, vjust=1, color="red")+theme_bw()

rda_ind
## Call: rda(formula = snp_pcs ~ as.matrix(env[, c(2)]) +
## as.matrix(env[, c(3)]) + as.matrix(env[, 4]) +
## (as.matrix(dbmem_pos[, 1])) + (as.matrix(dbmem_pos[, 2])) +
## (as.matrix(dbmem_pos[, 3])) + (as.matrix(dbmem_pos[, 4])) +
## (as.matrix(dbmem_pos[, 5])) + (as.matrix(dbmem_pos[, 6])) +
## (as.matrix(dbmem_pos[, 7])) + (as.matrix(dbmem_pos[, 8])))
## 
##               Inertia Proportion Rank
## Total          6.4115     1.0000     
## Constrained    4.9251     0.7682   11
## Unconstrained  1.4864     0.2318  103
## Inertia is variance 
## 
## Eigenvalues for constrained axes:
##  RDA1  RDA2  RDA3  RDA4  RDA5  RDA6  RDA7  RDA8  RDA9 RDA10 RDA11 
## 3.758 1.076 0.056 0.016 0.011 0.006 0.002 0.000 0.000 0.000 0.000 
## 
## Eigenvalues for unconstrained axes:
##     PC1     PC2     PC3     PC4     PC5     PC6     PC7     PC8 
## 0.11497 0.05451 0.02401 0.02194 0.02174 0.02118 0.02082 0.02073 
## (Showing 8 of 103 unconstrained eigenvalues)
aov_global
rda_ind_con_spatial
## Call: rda(formula = snp_pcs ~ as.matrix(env[, c(2)]) +
## Condition((as.matrix(dbmem_pos)[, c(1:8)])))
## 
##                 Inertia Proportion Rank
## Total         6.4115155  1.0000000     
## Conditional   4.8953346  0.7635222    8
## Constrained   0.0043957  0.0006856    1
## Unconstrained 1.5117852  0.2357922  103
## Inertia is variance 
## 
## Eigenvalues for constrained axes:
##     RDA1 
## 0.004396 
## 
## Eigenvalues for unconstrained axes:
##     PC1     PC2     PC3     PC4     PC5     PC6     PC7     PC8 
## 0.13284 0.05644 0.02721 0.02201 0.02176 0.02119 0.02086 0.02074 
## (Showing 8 of 103 unconstrained eigenvalues)
aov_con_spatial 
rda_ind_con_temp
## Call: rda(formula = snp_pcs ~ Condition(as.matrix(env[, c(2)])) +
## (as.matrix(dbmem_pos)[, c(1:8)]))
## 
##               Inertia Proportion Rank
## Total          6.4115     1.0000     
## Conditional    0.1808     0.0282    1
## Constrained    4.7189     0.7360    8
## Unconstrained  1.5118     0.2358  103
## Inertia is variance 
## 
## Eigenvalues for constrained axes:
##  RDA1  RDA2  RDA3  RDA4  RDA5  RDA6  RDA7  RDA8 
## 3.582 1.056 0.054 0.014 0.009 0.004 0.001 0.000 
## 
## Eigenvalues for unconstrained axes:
##     PC1     PC2     PC3     PC4     PC5     PC6     PC7     PC8 
## 0.13284 0.05644 0.02721 0.02201 0.02176 0.02119 0.02086 0.02074 
## (Showing 8 of 103 unconstrained eigenvalues)
aov_con_temp

notes for going forward
* it might be innappropirate to conduct anovas of rda condition on the combined environmental or dbmem matrices, it may be neccessary to indivudally test each variable, consitioned on all other variables… read into this before proceeding
* if we conduct the RDA on all pcs (not just those exceeding the kaiser guttman criteria) we can relate the RDA results more clearly to FST, i.e. we can see how much of the total variation among the cline is explained by RDA and compare this to results from amova and so on